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Abstract: 

<p>Time-reversal  is  a  method  of  focusing  sound  in  the  ocean  that  has  found  a  variety  of 
applications  in  recent  years,  ranging  from  underwater  communications  to  biological  stone 
destruction.  In  order  to  produce  a  focal  spot,  the  time-reversal  process  first  needs  to  acquire 
the  Green's  function  between  source  and  receiver.  This  Green's  function  is  time-reversed  and 
retransmitted  in  order  to  produce  a  spatio-temporal  focal  spot  at  the  original  source  location.  If 
the  medium  properties  in  between  the  source  and  receiver  change  between  the  acquisition  of  the 
Green's  function  and  the  subsequent  retransmission,  the  quality  of  the  focal  spot  can  degrade 
or  even  disappear.  However,  the  time-reversal  focal  spot  has  been  found  to  be  surprisingly 
robust  to  changes  in  medium  properties,  which  are  chiefly  sound  speed  fluctuations  in  underwater 
acoustics.  At  445  Hz,  the  focal  spot  was  seen  to  persist  for  a  week,  while  at  3.5  kHz,  the  focal  spot 
persists  for  about  an  hour.</p>  <p>Sensitivity  kernels  have  the  ability  to  linearly  map  sound  speed 
perturbations  to  a  perturbation  of  an  acoustic  parameter  such  as  travel-time  or  pressure.  Sensitivity 
kernels  have  a  Fresnel-like  interference  pattern  with  regions  of  positive  and  negative  sensitivities 
in  the  medium.  Time-reversal,  which  causes  different  arrival  paths  to  arrive  at  the  same  time, 
results  in  overlapping  sensitivity  kernels  that  leads  to  a  net  reduction  in  pressure  sensitivity  at 
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ABSTRACT  OF  THE  DISSERTATION 


Pressure  Sensitivity  Kernels  Applied  to  Time-reversal  Acoustics 

by 

Kaustubha  Raghukumar 
Doctor  of  Philosophy  in  Oceanography 

University  of  California  San  Diego,  2009 
William  S.  Hodgkiss,  Chair 

Time-reversal  is  a  method  of  focusing  sound  in  the  ocean  that  has  found  a 
variety  of  applications  in  recent  years,  ranging  from  underwater  communications 
to  biological  stone  destruction.  In  order  to  produce  a  focal  spot,  the  time- reversal 
process  hrst  needs  to  acquire  the  Green’s  function  between  source  and  receiver. 
This  Green’s  function  is  time-reversed  and  retransmitted  in  order  to  produce  a 
spatio-temporal  focal  spot  at  the  original  source  location.  If  the  medium  prop¬ 
erties  in  between  the  source  and  receiver  change  between  the  acquisition  of  the 
Green’s  function  and  the  subsequent  retransmission,  the  quality  of  the  focal  spot 
can  degrade  or  even  disappear.  However,  the  time-reversal  focal  spot  has  been 
found  to  be  surprisingly  robust  to  changes  in  medium  properties,  which  are  chiefly 
sound  speed  fluctuations  in  underwater  acoustics.  At  445  Hz,  the  focal  spot  was 
seen  to  persist  for  a  week,  while  at  3.5  kHz,  the  focal  spot  persists  for  about  an 
hour. 

Sensitivity  kernels  have  the  ability  to  linearly  map  sound  speed  perturba¬ 
tions  to  a  perturbation  of  an  acoustic  parameter  such  as  travel-time  or  pressure. 
Sensitivity  kernels  have  a  Fresnel-like  interference  pattern  with  regions  of  positive 
and  negative  sensitivities  in  the  medium.  Time-reversal,  which  causes  different  ar¬ 
rival  paths  to  arrive  at  the  same  time,  results  in  overlapping  sensitivity  kernels  that 
leads  to  a  net  reduction  in  pressure  sensitivity  at  the  focal  spot.  Upon  expressing 
the  pressure  at  the  focal  spot  in  terms  of  sensitivity  kernels,  source  transmissions 
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are  derived  that  are  even  more  robust  than  time-reversal. 

The  theory  developed  using  pressure  sensitivity  kernels  is  tested  on  experi¬ 
mental  data,  along  with  an  internal  wave  model,  using  various  metrics.  The  linear 
limitations  of  the  kernels  are  explored  in  the  context  of  time-evolving  Green’s  func¬ 
tions.  The  optimized  source  functions  are  then  tested  using  experimental  Green’s 
functions  and  their  behavior  is  seen  to  be  in  the  right  sense. 

Finally,  thermistor  chain  data  gathered  during  a  similar  experiment  allowed 
for  the  testing  of  a  ’synthetic  aperture  thermistor  chain’.  Temperature  observa¬ 
tions,  presumably  caused  by  displacements  of  a  reference  temperature  prohle,  are 
inverted  to  provide  an  estimate  of  the  background  temperature  profile. 


Chapter  1 


Introduction 


1.1  Background 

The  scientific  focus  on  shallow  water  acoustics  since  the  end  of  the  cold  war 
has  lead  to  a  renewed  interest  in  environmentally-induced  acoustic  fluctuations, 
together  with  ways  and  means  to  understand  and  mitigate  their  effect.  A  wide 
variety  of  applications  ranging  from  underwater  communications  and  mine  detec¬ 
tion  to  the  acoustic  tagging  of  marine  mammals  require  an  understanding  of  the 
oceanographic  fluctuations  that  can  influence  an  acoustic  signal  as  it  makes  its  way 
from  a  source  to  the  receiver.  Shallow  water  acoustic  propagation,  in  particular, 
has  to  deal  with  a  higher  degree  of  environmental  complexity  compared  to  acoustic 
propagation  in  deeper  water.  This  is  due  to  a  variety  of  reasons  such  as  increased 
bottom  and  surface  interaction,  a  higher  frequency  of  operation,  and  a  higher  noise 
floor  due  to  the  proximity  to  shipping  lanes.  [1]  Another  major  cause  of  fluctua¬ 
tions  are  internal  waves  whose  presence  results  in  the  heaving  of  the  thermocline 
in  the  water  column,  resulting  in  sound  speed  perturbations  that  in  turn  perturb 
the  acoustic  signal. 

Understanding  the  effects  of  a  fluctuating  sound  speed  prohle  on  the  var¬ 
ious  paths  of  the  received  multipath  signal  is  important  in  order  to  separate  out 
environmental  effects  from  the  desired  signal.  In  order  to  be  able  to  better  perform 
physics-based  signal  processing,  a  variety  of  experiments  have  been  conducted  in 
recent  years  to  study  acoustic  fluctuations.  As  an  example,  one  such  cruise  that 


1 


2 


took  wide-ranging  environmental  measnrements  was  the  Kani-Ex  experiment  that 
took  place  off  the  coast  of  Kanai  in  July  2003.  [2]  In  addition  to  near-continuous 
probing  of  the  acoustic  channel  using  a  variety  of  vertical  and  bottom  mounted 
arrays,  oceanographic  measurements  were  made  of  quantities  such  as  the  surface 
wave  spectrum,  2D  temperature  structure,  salinity,  currents  and  bottom  proper¬ 
ties. 

Complementary  to  the  work  done  during  Kauai-Ex,  which  can  primarily  be 
seen  as  dealing  with  understanding  the  forward  problem,  is  the  effort  to  invert  an 
acoustic  signal  for  the  environment  that  gave  it  its  unique  signature.  This  work 
ranges  from  the  contributions  to  geoacoustic  inversion  made  by  Gerstoft  [3]  to  the 
shallow  water  acoustic  tomography  work  by  Cornuelle  et.  ah  [4]  and  Hursky  et.  ah 

[5] 

Time-reversal  is  a  method  of  focusing  sound  in  the  ocean  in  range,  depth 
and  time  that  has  proven  remarkably  robust  to  sound  speed  fluctuations.  While 
time-reversal  has  been  widely  used  in  applications  ranging  from  underwater  com¬ 
munications  to  target  detection,  efforts  to  understand  the  reasons  for  its  stability 
have  been  lacking.  Fundamentally,  this  thesis  can  be  seen  as  providing  an  explana¬ 
tion  for  the  stability  of  time-reversal  using  some  of  the  tools  developed  in  seismic 
and  acoustic  tomography.  Further,  the  same  tools  are  utilized  in  an  optimization 
framework  together  with  a  knowledge  of  fluctuation  statistics  to  develop  new  sets 
of  signals  that  can  be  made  even  more  stable  than  time-reversal.  This  can  be  seen 
as  an  example  of  physics-based  signal  processing  that  was  alluded  to  earlier. 

1.2  Time-reversal 

The  hrst  time-reversal  experiment  in  the  ocean  was  conducted  by  Parvulescu 
and  Clay  in  1965.  [6,  7]  The  authors  viewed  the  ocean  as  a  black  box  that  has  a 
particular  impulse  response.  Given  an  input  signal  s{t)  and  an  impulse  response 
h(t),  the  received  signal  r(t)  can  be  expressed  as  the  convolution  of  the  input  signal 
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with  the  impulse  response. 


r{t) 


s{t)  ®  hit) 


s{T)h{t 


T)dT 


(1.1) 

(1.2) 


Recognizing  that  a  convolution  operation  with  one  of  the  signals  reversed  in  time 
is  the  same  as  correlating  the  two  signals,  they  were  able  to  use  the  ocean  as  a 
‘correlator’  by  time-reversing  and  retransmitting  an  approximation  to  the  impulse 
response,  i.e.  by  recording  and  transmitting  h{—t),  they  were  able  to  obtain  a 
received  signal  that  is  the  autocorrelation  of  h(t).  Of  course,  there  is  the  possibility 
that  the  impulse  response  may  have  changed  between  the  time  of  recording  and 
retransmitting,  in  which  case  the  received  signal  will  be  the  cross-correlation  of 
two  impulse  responses.  The  correlation  peak  is  akin  to  a  focal  spot  at  the  location 
of  the  sonrce,  whose  qnality  (peak  power,  peak-to-sidelobe  ratio)  will  depend  on 
the  degree  of  similarity  between  the  two  impulse  responses.  Fig.  1.1,  taken  from 
Parvnlescn[7],  is  a  simple  schematic  of  how  time-reversal  works.  Time-reversal  of 
the  impnlse  response  and  snbseqnent  retransmission  will  cause  each  pulse  to  elicit 
a  response.  One  of  the  paths  for  each  response  will  arrive  at  exactly  the  same  time, 
and  dne  to  the  linearity  of  the  hydrophone  response,  the  pressure  at  this  time  will 
be  exactly  three  times  those  of  the  sidelobes.  One  can  thus  see  that  the  larger  the 
number  of  multipaths,  the  higher  will  be  the  peak  relative  to  the  sidelobes. 

Time-reversal  is  a  time-domain  concept  whose  eqnivalent  in  the  freqnency 
domain  is  phase  conjugation.  Phase  conjngation  was  originally  employed  in  non¬ 
linear  optics  where  medinm  nonlinearities  are  used  to  scatter  incident  light  in  a 
way  that  the  resulting  scattered  signal  is  a  phase-conjngated  version  of  the  incident 
radiation.  [8]The  scattering  medinm  is  called  a  phase  conjngating  mirror,  and  jnst 
like  in  aconstics,  takes  advantage  of  the  time-reversal  invariance  of  the  wave  eqna- 
tion.  In  addition  to  the  invariance  property,  time-reversal  also  takes  advantage  of 
reciprocity  in  the  medinm,  i.e.  the  response  of  the  medinm  is  the  same  whether 
going  from  sonrce  to  receiver  or  receiver  to  sonrce. 

Following  a  nnmber  of  years  of  attention  in  the  optics  commnnity,  phase 
conjngation  regained  prominence  in  acoustics  after  a  series  of  ultrasonics  exper¬ 
iments  by  Fink  et.  ah  [9,  10]  who  demonstrated  focusing  of  sonnd  using  time 


Receive 


4 


Transmit 

impulse 

— - ►  t 


Time-reversed 

signal 

h* —  T  — H 


1 


2 


3 


♦  t 


Travel  path 
1 


2 


3 


Total  signal 


peak 

-2T- 


Figure  1.1:  A  schematic  for  how  time-reversal  works,  based  on  Parvulescu  and 
Clay.  [6] 

reversal  mirrors.  In  one  of  the  hrst  experiments,  a  closed  cavity  containing  a 
sonrce  is  snrronnded  by  transdncers.  Sonnd  emitted  by  the  sonrce  is  scattered  by 
heterogeneities  in  the  cavity  and  recorded  by  the  snrronnding  transdncers.  Upon 
emitting  the  time-reversed  version  of  the  recording,  it  was  seen  to  focns  back  on 
the  sonrce  location.  This  experiment  showed  promise  in  the  ability  of  sonnd  fo¬ 
cused  by  time-reversal  in  biological  stone  destruction.  Another  experiment  showed 
that  if  the  medium  is  weakly  homogeneous  (dominated  by  single  scattering),  then  a 
time-reversal  mirror  can  produce  a  focal  spot  comparable  to  that  of  a  closed- cavity, 
whose  realization  can  be  difficult  in  practice. 


5 


Around  the  same  time  that  Fink  et.  ah  were  conducting  the  ultrasonic 
time-reversal  experiments,  Jackson  and  Dowling  [11]  proposed  the  viability  of 
phase  conjugation  as  a  means  of  solving  many  longstanding  problems  in  under¬ 
water  acoustics.  The  behavior  of  phase  conjugate  arrays  in  underwater  acoustics 
was  studied  through  several  examples  and  the  effects  of  aperture  size  and  inho¬ 
mogeneities  in  medium  propagation  were  treated  for  both  near-held  and  far-held 
regions. 

The  hrst  time-reversal  array  experiment  in  the  ocean  was  conducted  by 
Kuperman  et.  al.  [12]  in  the  Mediterranean  Sea  in  1996.  This  experiment  suc¬ 
cessfully  demonstrated  focusing  in  the  ocean  waveguide  at  a  distance  of  6.3  km 
using  a  50  ms  pulse  centered  at  445  Hz.  In  addition  to  showing  that  the  ocean 
acoustic  experiments  were  consistent  with  theory,  the  stability  of  the  focal  spot 
was  studied,  and  found  to  be  robust  for  a  time  period  of  over  an  hour.  The  success 
of  this  proof-of-concept  experiment  in  1996  paved  the  way  for  a  series  of  time- 
reversal  experiments  in  the  region  conducted  by  the  same  group,  and  were  termed 
the  Focused  Acoustic  Fields  (FAF)  experiments.  These  experiments  [13,  14,  15,  16] 
studied  the  utility  of  time-reversal  in  a  variety  of  applications,  ranging  from  un¬ 
derwater  communications  to  buried  target  detection.  Some  of  these  applications 
will  be  described  in  the  following  section.  In  addition,  the  experiments  also  exam¬ 
ined  the  stability  of  the  focal  spot  for  a  number  of  ranges,  frequencies  and  depths, 
providing  valuable  insight  into  the  physics  of  time-reversal  in  the  ocean.  The  data 
gathered  during  one  of  these  experiments,  FAF  2005,  is  used  in  this  thesis  while 
studying  the  stability  aspects  of  the  time-reversal  focus. 

1.2.1  Applications  of  Time-reversal 

The  ability  of  time-reversal  to  focus  sound  in  space  and  time  lends  its  itself 
to  a  wide  variety  of  applications.  Time-reversal,  when  proposed  by  Parvulescu  and 
Clay  [6],  was  seen  as  a  way  of  allowing  the  ocean  to  perform  the  correlation  process, 
thus  saving  on  the  limited  computational  resources  available  at  the  time.  When 
phase  conjugation  was  later  proposed  in  optics,  it  was  seen  as  having  potential 
applications  in  adaptive  optics,  lenseless  imaging,  phase-conjugate  resonators  and 
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image  processing[17]. 

Time-reversal  in  ultrasonics [10]  was  originally  proposed  as  a  way  of  per¬ 
forming  lithotripsy  (destruction  of  biological  stones),  by  virtue  of  being  able  to 
focus  sound  inside  a  cavity.  In  addition,  it  was  found  that  iterative  time-reversal 
in  the  presence  of  multiple  scatterers  tends  to  converge  to  a  focus  on  the  strongest 
scatterer.  This  aspect  of  time-reversal,  while  also  having  biomedical  applications, 
has  been  used  in  non-destructive  testing  to  detect  small  surface  anomalies. 

In  underwater  acoustics,  the  most  popular  use  of  time-reversal  has  been  in 
the  area  of  underwater  communications.  The  idea  was  hrst  proposed  by  Roux  [18] , 
who  recognized  that  time-reversal  has  the  ability  to  mitigate  the  effects  of  multi- 
path,  which  have  long  plagued  underwater  communications.  In  communications, 
where  several  symbols  are  transmitted  one  after  the  other,  multipath  causes  the 
‘echos’  associated  with  any  one  symbol  to  overlap  with  the  other  symbols,  causing 
what  is  termed  inter-symbol  interference  (ISI).  Typically,  communication  receivers 
are  built  to  remove  ISI  using  signal  processing  tools  known  as  equalizers [19].  Time- 
reversal,  on  the  other  hand  exploits  multipath,  by  effectively  delaying  each  multi- 
path  such  that  all  the  different  paths  arrive  at  the  same  time,  causing  a  focal  spot. 
Edelmann[14]  showed  that  it  is  possible  to  encode  phase  information  in  the  time- 
reversed  Green’s  function  that  can  then  be  recovered  at  the  receiver  end.  The  use 
of  the  time-reversed  Green’s  function  as  a  communication  symbol  has  the  advan¬ 
tage  of  signihcantly  reduced  ISI  due  to  the  creation  of  a  focal  spot  for  each  symbol. 
The  ISI  suppression  achieved  by  time-reversal  is  often  insufficient  for  higher  data 
rate  communications,  or  multiple-input  multiple-output  (MIMO)  communications 
[20],  in  which  case  time-reversal  is  typically  followed  by  receiver  side  equalization. 

Other  uses  of  time-reversal  in  underwater  acoustics  include  its  use  as  an 
acoustic  barrier  [21]  and  to  detect  buried  objects[22].  In  the  acoustic  barrier  ex¬ 
periment,  a  disturbance  in  the  focal  plane  is  used  to  alert  the  presence  of  a  foreign 
object  between  the  source  and  receiver.  When  used  for  target  detection,  the  ability 
of  time-reversal  to  focus  on  a  target  results  in  less  energy  interacting  with  rough 
ocean  boundaries  above  and  below  the  focal  spot.  This  give  rise  to  an  enhanced 
echo-to-reverberation  ratio  which  can  be  exploited  to  localize  the  target. 
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1.2.2  Stability  of  Time-reversal 

A  highly  attractive  aspect  of  time-reversal,  besides  its  focusing  ability,  is 
the  fact  that  the  focus  produced  by  a  particular  Green’s  function  can  persist  for 
a  while.  This  happens  despite  the  fact  that  the  Green’s  function  between  the 
acquisition  and  subsequent  retransmissions  has  changed.  At  lower  frequencies  (445 
Hz),  it  was  found  that  the  focal  spot  persists  for  several  hours  upto  a  week  [13].  At 
higher  frequencies  (3.5  kHz),  the  focal  spot  persisted  for  about  an  hour  [15].  The 
authors  examined  fluctuations  of  the  time-reversal  focal  spot  using  the  waveguide 
invariant.  They  found  that,  in  the  presence  of  internal  wave  induced  fluctuations, 
the  focal  spot  merely  shifted  in  range,  but  was  still  present  at  a  different  location. 
While  the  explanation  is  valid  for  a  single  frequency,  or  the  narrowband  case,  the 
translation  to  a  broadband  scenario  is  not  entirely  obvious.  Recent  experiments 
(FAF  2005)  were  carried  out  using  a  1  kHz  bandwidth  centered  around  3.5  kHz. 
Various  other  explanations  have  also  been  put  forth  to  explain  the  stability  of  the 
focus,  which  are  examined  in  more  detail  in  Ghapter  2. 

1.3  Motivation  and  purpose 

This  thesis  can  be  seen  as  an  attempt  to  explain  the  stability  of  the  time- 
reversal  focal  spot  that  can  be  applied  regardless  of  the  specific  characteristics 
of  the  time-reversal  experiment  such  as  geometry,  frequency,  bandwidth,  etc.  One 
constraint  that  is  invoked  is  that  the  sound  speed  perturbations  that  cause  a  degra¬ 
dation  of  the  focal  spot  be  small  enough  to  satisfy  the  single  scattering  assumption 
(first  Born  approximation).  It  is  also  assumed  that  changes  to  the  focal  spot  are 
caused  solely  by  sound  speed  fluctuations,  and  the  effects  of  features  such  as  sur¬ 
face  waves  and  non-linear  internal  waves  are  negligible.  This  is  because  non-linear 
internal  waves  have  been  observed [23]  to  cause  horizontal  refraction  of  acoustic  en¬ 
ergy.  While  the  analysis  in  this  thesis  can  easily  be  extended  to  three  dimensions, 
only  two-dimensional  effects  are  considered  at  this  time. 

A  change  in  the  pressure  of  the  focal  spot  due  to  a  change  in  the  sound  speed 
profile  can  be  thought  of  as  a  pressure  sensitivity  to  sound  speed  perturbations. 


Time-reversal  has  been  observed  to  have  a  lower  pressure  sensitivity  than  a  simple 
pulse  transmission.  Given  this  line  of  reasoning,  it  becomes  apparent  that  the  use 
of  a  formal  sensitivity  analysis  framework  can  yield  potentially  useful  results. 

The  FAF  experiment  that  was  conducted  in  July  2005  in  the  Mediterranean 
sea  focused  on  several  aspects  of  time-reversal  and  time-reversal  communications 
using  vertical  line  arrays  and  billboard  arrays.  The  focal  spot  was  seen  to  persist  for 
about  an  hour  at  3.5  kHz  with  a  1  kHz  pulse  bandwidth.  In  addition  to  the  time- 
reversal  segment,  the  experiment  also  contained  a  portion  where  Green’s  functions 
were  acquired  every  20  seconds  over  8  hours.  Given  the  availability  of  Green’s 
functions,  the  time-reversal  focus  can  be  passively  constructed  on  the  computer, 
and  a  variety  of  focusing  scenarios  can  be  studied. 

In  addition  to  the  rich  acoustic  dataset  available,  the  experiment  also  gath¬ 
ered  GTD  (conductivity-temperature-depth)  prohles  and  temperature  time-series 
over  12  days  using  two  moored  thermistor  chains.  The  availability  of  acoustic  and 
environmental  data  allows  for  a  detailed  study  of  fluctuations  and  their  effect  on 
time-reversal. 

The  main  questions  that  this  thesis  addresses  can  thus  be  summarized  as 

follows: 

•  Does  time-reversal  have  a  lower  pressure  sensitivity  to  sound  speed  fluctua¬ 
tions  than  a  one-way  transmission? 

•  Gan  sensitivity  analysis  be  used  to  provide  an  answer  why? 

•  Gan  knowledge  of  the  pressure  sensitivity  be  exploited  to  develop  signals  that 
are  even  less  sensitive  than  time-reversal? 

•  Gan  knowledge  of  the  pressure  sensitivity  be  used  to  predict  the  behavior  of 
the  time-reversal  focal  spot? 

The  availability  of  the  thermistor  chain  record  also  allowed  for  an  additional 
investigation  that  could  be  useful  in  sensitivity  studies  and  internal  wave  model¬ 
ing.  The  hypothesis  is  that  the  temperature  fluctuations  observed  are  caused  by 
repeated  displacements  of  a  background,  or  reference  temperature  prohle.  Thus, 
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given  the  thermistor  chain  record,  it  should  be  possible  to  invert  for  the  reference 
prohle  using  iterative  methods.  The  last  question  question  addressed  in  this  thesis 
is  therefore, 

•  Given  a  thermistor  chain  record,  is  it  possible  to  invert  for  a  background 
reference  temperature  prohle? 


1.4  Outline 

This  thesis  is  organized  as  follows^.  Chapter  2  discusses  the  theory  for  sensi¬ 
tivity  analysis  applied  to  time-reversal  acoustics  by  introducing  pressure  sensitivity 
kernels.  [24]  Two  simple  cases  are  considered:  free-space  and  a  Pekeris  waveguide. 
Sensitivity  kernels  are  illustrated  for  each  case,  and  an  argument  is  made  for  the 
reduction  of  sensitivity  when  time-reversal  is  used.  The  theory  for  pressure  sensitiv¬ 
ity  kernels  is  then  taken  a  step  further,  where  the  knowledge  of  pressure  sensitivity 
is  used  to  derive  optimized  source  functions  that  are  closely  related  to  the  inverse 
hlter.[25]  These  functions  contain  adjustable  parameters  that  allow  for  a  tradeoff 
between  focal  spot  quality  and  sensitivity  to  sound  speed  huctuations.  Various 
features  of  these  functions  are  then  illustrated  for  three  simple  perturbation  sce¬ 
narios:  hrst,  when  medium  perturbations  are  uncorrelated,  second  when  there  is  a 
channel  of  zero  perturbation  variance  in  depth,  and  the  third  when  perturbations 
have  a  long  decorrelation  scale. 

Chapter  3  takes  the  theory  developed  in  Chapter  2  and  applies  it  to  the 
FAF05  dataset.  [26]  An  internal  wave  model  is  developed  whose  characteristics 
are  matched  to  the  environmental  data  available.  It  was  found  that  the  internal 
wave  model,  together  with  a  parabolic  equation  acoustic  code,  was  able  to  capture 
the  of  evolution  of  the  peak  of  the  time-reversal  focal  spot.  This  can  be  seen  as 
a  acoustic  test  of  the  internal  wave  model.  Sensitivity  kernels  are  then  computed 
using  experimentally  gathered  sound  speed  prohles  and  an  attempt  is  made  to  track 
Green’s  functions  as  the  evolve  using  these  kernels.  However,  it  was  found  that  the 

^Each  chapter  is  essentially  a  full  or  partial  reprint  of  papers  that  have  been  published, 
accepted,  submitted  or  in  preparation  for  submission  to  a  professional  journal.  Therefore,  each 
has  been  written  in  a  paper  format  and  is  self-contained. 
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linear  nature  of  the  kernels  limits  the  tracking  ability  to  only  about  10  minutes  of 
fluctuations  for  the  perturbation  scenario  encountered.  Optimized  source  functions 
are  then  computed  using  the  sensitivity  kernels  and  the  experimentally  gathered 
Green’s  functions.  Their  performance  is  compared  to  time-reversal  and  the  inverse 
hlter  for  the  given  dataset  and  found  to  be  in  the  right  sense. 

Chapter  4  deals  with  the  problem  of  temperature  profile  estimation  using 
the  thermistor  chain  record.  Thermistor  chain  data  used  belonged  to  the  FAF06 
experiment.  A  different  set  of  thermistor  chain  data  is  used  compared  to  Chapter 
3  due  to  the  closer  depth  and  time-spacing  available  in  the  FAF06  thermistor 
chain  data.  This  allowed  for  a  better  demonstration  of  the  estimation  method. 
The  iterative  Gauss-Newton  method  that  is  used  is  introduced  and  its  application 
to  the  problem  is  discussed.  The  prohle  estimation  algorithm  is  then  applied  to 
three  different  cases.  The  hrst  two  use  synthetically  generated  data  that  have  a 
known  reference  temperature  while  the  third  uses  the  FAF06  temperature  data  to 
estimate  an  unknown  reference.  The  first  synthetic  data  case  study  uses  randomly 
generated  temperature  samples,  while  the  second  uses  the  internal  wave  model  to 
generate  temperature  time-series. 


Bibliography 


[1]  W.  A.  Kuperman  and  J.  F.  Lynch.  Shallow-water  acoustics.  Physics  Today, 
57:55-61,  2004. 

[2]  M.  B.  Porter,  P.  Hursky,  M.  Siderius,  M.  Badiey,  J.  Caruthers,  W.  S.  Hodgkiss, 
K.  Raghukumar,  D.  Rouseff,  W.  Fox,  C.  Moustier,  B.  Calder,  B.  J.  Kraft, 
K.  McDonald,  P.  Stein,  J.  K.  Lewis,  and  S.  Rajan.  Kauaiex:  Environmen¬ 
tal  effects  on  hf  acoustic  propagation  with  application  to  communications. 
J.  Acoustic  Soc.  Am.,  115:2468-2468,  2004. 

[3]  P.  Gerstoft.  Inversion  of  seismo-acoustic  data  using  genetic  algorithms  and  a 
posteriori  probability  distributions,.  J.  Acoustic  Soc.  Am.,  95:770-782,  1994. 

[4]  B.  Cornuelle,  M.  Dzieciuch,  and  E.  Skarsoulis.  Wave-theoretic  kernels  for 
long-range  ocean  travel-time  measurement.  J.  Acoustic  Soc.  Am.,  119:3248, 
2006. 

[5]  P.  Hursky,  M.  B.  Porter,  B.  D.  Cornuelle,  W.  S.  Hodgkiss,  and  W.  A.  Ku¬ 
perman.  Adjoint  modelling  for  acoustic  inversion.  J.  Acoustic  Soc.  Am., 
115:607-619,  2004. 

[6]  A.  Parvulescu  and  C.  S.  Clay.  Reproducibility  of  signal  transmissions  in  the 
ocean.  Radio  Electron  Eng.,  29:223-228,  1965. 

[7]  A.  Parvulescu.  Matched-signal  (“MESS”)  processing  by  the  ocean.  J.  Acous- 
tic  Soc.  Am.,  98(2):943-960,  1995. 

[8]  B.  Y.  ZehDovich,  N.  F.  Pilipetsky,  and  V.  V.  Shkunov.  Principles  of  phase 
conjugation.  Springer- Verlag,  Berlin,  1985. 

[9]  M  Fink,  C.  Prada,  F.  Wu,  and  D.  Cassereau.  Self-focusing  with  time  reversal 
mirror  in  inhomogeneous  media.  Proc.  IEEE  Ultrason.  Symp.  1989  Montreal, 
2:681-686, 1989. 

[10]  M.  Fink.  Time-reversal  mirrors.  Journal  of  Physics  D:  Applied  Physics, 
26(9):1333-1350,  1993. 


11 


12 


[11]  D.  R.  Jackson  and  D.  R.  Dowling.  Phase  conjugation  in  underwater  acoustics. 
J.  Acoustic  Soc.  Am.,  89:171-181,  1991. 

[12]  W.  A.  Kuperman,  W.  S.  Hodgkiss,  H.  C.  Song,  T.  Akal,  C.  Ferla,  and  D.  R. 
Jackson.  Phase  conjugation  in  the  ocean:  Experimental  demonstration  of  an 
acoustic  time-reversal  mirror.  J.  Acoustic  Soc.  Am.,  103:25-40,  1998. 

[13]  W.  S.  Hodgkiss,  H.  C.  Song,  W.  A.  Kuperman,  T.  Akal,  C.  Ferla,  and  D.  R. 
Jackson.  A  long-range  and  variable  focus  phase-conjugation  experiment  in 
shallow  water.  J.  Acoustic  Soc.  Am.,  105:1597-1604,  1999. 

[14]  G.  F.  Edelmann,  T.  Akal,  W.  S.  Hodgkiss,  S.  Kim,  W.  A.  Kuperman,  and 
H.  C.  Song.  An  initial  demonstration  of  underwater  acoustic  communication 
using  time-reversal.  IEEE  J.  Oceanic  Engineering,  27:602-608,  2002. 

[15]  S.  Kim,  W.  A.  Kuperman,  W.  S.  Hodgkiss,  H.  C.  Song,  and  G.  F.  Edelmann. 
Robust  time-reversal  focusing  in  the  ocean.  J.  Acoustic  Soc.  Am.,  114:145- 
157,  2001. 

[16]  H.  G.  Song,  W.  A.  Kuperman,  W.  S.  Hodgkiss,  T.  Akal,  and  G.  Ferla.  Iterative 
time  reversal  in  the  ocean.  J.  Acoustic  Soc.  Am.,  105:3176-3184,  1999. 

[17]  D.  M.  Pepper.  Applications  of  optical  phase  conjugation.  Scientific  American, 
254:74-83,  1986. 

[18]  P.  Roux.  Time-reversal  mirrors:  Application  to  waveguide  propagation  and 
characterization  of  fluid  flows  (French).  Phd  thesis,  University  of  Paris,  1997. 

[19]  M.  Stojanovic.  Recent  advances  in  high-speed  underwater  acoustic  communi¬ 
cations.  IEEE  J.  Oceanic  Engineering,  21:125-136,  1996. 

[20]  H.  G.  Song,  W.  S.  Hodgkiss,  W.  A.  Kuperman,  T.  Akal,  and  M.  Steven¬ 
son.  Multiuser  communications  using  passive  time  reversal.  IEEE  J.  Oceanic 
Engineering,  32:915-926,  2007. 

[21]  H.  G.  Song,  W.  A.  Kuperman,  W.  S.  Hodgkiss,  T.  Akal,  and  P.  Guerrini. 
Demonstration  of  a  high-frequency  acoustic  barrier  with  a  time-reversal  mir¬ 
ror.  IEEE  J.  Oceanic  Engineering,  28:246-249,  2003. 

[22]  S.  Kim,  W.  A.  Kuperman,  W.  S.  Hodgkiss,  H.  G.  Song,  G.  F.  Edelmann,  and 
T.  Akal.  Echo-to-reverberation  enhancement  using  a  time-reversal  mirror. 
J.  Acoustic  Soc.  Am.,  115(4):1525-1531,  2004. 

[23]  S.  D.  Frank,  M.  Badiey,  J.  F.  Lynch,  and  W.  L.  Siegmann.  Experimental  evi¬ 
dence  of  three-dimensional  acoustic  propagation  caused  by  nonlinear  internal 
waves.  JASA,  118:723-734,  2005. 


13 


[24]  K.  Raghukumar,  Bruce  D.  Cornuelle,  William  S.  Hodgkiss,  and  William  A. 
Kuperman.  Pressure  sensitivity  kernels  applied  to  time-reversal  acoustics. 
J.  Acoustic  Soc.  Am.,  124:98-112,  2008. 

[25]  M.  Tauter,  J.  L.  Thomas,  and  M.  Fink.  Time-reversal  and  the  inverse  hlter. 
J.  Acoustic  Soc.  Am.,  108:223-234,  2000. 

[26]  K.  Raghukumar,  Bruce  D.  Cornuelle,  William  S.  Hodgkiss,  and  William  A. 
Kuperman.  Experimental  demonstration  of  pressure  sensitivity  kernels. 
J.  Acoustic  Soc.  Am.,  2009,  to  be  submitted. 


Chapter  2 


Pressure  Sensitivity  Kernels 
Applied  to  Time-reversal 
Acoustics 


[Reprinted  with  permission  from  Kaustubha  Raghukumar,  Bruce  Cornuelle, 
William  Hodgkiss,  William  Kuperman,  J.  Acoustic  Soc.  Am.,  Vol.  124,  Issue  1, 
pp.  98-112,  2008,  Acoustical  Society  of  America] 

Sensitivity  kernels  for  receptions  of  broadband  sound  transmissions  are  used 
to  study  the  effect  of  the  transmitted  signal  on  the  sensitivity  of  the  reception  to 
environmental  perturbations.  A  first-order  Born  approximation  is  used  to  obtain 
the  pressure  sensitivity  of  the  received  signal  to  small  changes  in  medium  sound 
speed.  The  pressure  perturbation  to  the  received  signal  caused  by  medium  sound 
speed  changes  is  expressed  as  a  linear  combination  of  single-frequency  sensitivity 
kernels  weighted  by  the  signal  in  the  frequency  domain.  This  formulation  can  be 
used  to  predict  the  response  of  a  source  transmission  to  sound  speed  perturba¬ 
tions.  The  stability  of  time-reversal  is  studied  and  compared  to  that  of  a  one-way 
transmission  using  sensitivity  kernels.  In  the  absence  of  multipath,  a  reduction 
in  pressure  sensitivity  using  time-reversal  is  only  obtained  with  multiple  sources. 
This  can  be  attributed  both  to  the  presence  of  independent  paths  and  to  cancella¬ 
tions  that  occur  due  to  the  overlap  of  sensitivity  kernels  for  different  source-receiver 
paths.  The  sensitivity  kernel  is  then  optimized  to  give  a  new  source  transmission 


14 


15 


scheme  that  takes  into  account  knowledge  of  the  medium  statistics  and  is  related 
to  the  regularized  inverse  hlter. 

2.1  Introduction 

Time-reversed  acoustics,  due  to  its  ability  to  focus  energy  in  space  and  time, 
has  found  numerous  applications  ranging  from  underwater  and  wireless  commu¬ 
nications  to  medical  ultrasound  therapy.  [1,  2,  3]  Time-reversal  works  by  trans¬ 
mitting  a  signal  which  is  a  time-reversed  version  of  the  observed  Green’s  function 
between  a  source  and  receiver.  This  scheme  assumes  that  the  Green’s  function  has 
not  signihcantly  changed  between  its  acquisition  and  time-reversed  retransmission. 
Despite  this  assumption,  the  ability  of  time-reversal  to  produce  a  focal  spot  has 
been  found  to  be  remarkably  robust  to  environmental  changes.  In  the  ocean,  the 
focal  spot  has  been  shown  to  persist  for  a  week  at  445  Hz  (50  ms  pulse)  [4],  and 
close  to  an  hour  at  3.5  kHz  (10  ms  pulse).  [5] 

Various  explanations  have  been  offered  to  understand  the  reason  for  the 
robustness  of  time- reversal  to  environmental  perturbations.  Kim  et  al.  [5],  using  the 
theory  of  waveguide  invariants  for  an  oceanic  waveguide,  showed  that  while  there  is 
little  degradation  in  the  intensity  of  the  focal  spot,  it  is  seen  to  shift  back  and  forth 
in  range  as  a  result  of  internal  wave  induced  sound  speed  fluctuations.  While  any 
one  realization  of  the  focal  spot  may  or  may  not  be  stable,  an  ensemble  average  of 
focal  spots  tends  to  have  high  stability.  [6]  The  presence  of  spatial  diversity  (large 
arrays)  has  the  same  effect  as  an  ensemble  average  and  results  in  stability  of  the 
focal  spot.  In  random  media[7],  the  probability  of  a  realization  of  the  focal  spot 
differing  signihcantly  from  its  expected  value  is  small.  This  feature  was  referred 
to  as  the  self-averaging  property  of  the  focal  spot  which  leads  to  the  statistical 
stability  of  time-reversal  for  broadband  signals.  The  self-averaging  property  is 
based  on  specihc  asymptotic  limits  more  suited  to  optical  and  infrared  time-reversal 
than  sonar  or  ultrasound. 

The  above  contributions  suggest  the  use  of  a  formal  sensitivity  analysis  to 
analyze  the  stability  of  the  time-reversal  focal  spot  with  respect  to  sound  speed 
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perturbations.  Sensitivity,  when  defined  as  the  derivative  of  the  received  pressure 
signal  with  respect  to  medium  sound  speed  fluctuations,  provides  a  way  to  compare 
time-reversal  transmissions  to  one-way  transmissions.  Sensitivity  kernels  can  be 
a  useful  tool  for  analyzing  various  transmission  schemes  that  have  been  proposed 
that  seek  to  provide  an  even  more  stable  focal  spot  [5],  or  a  focal  spot  with  reduced 
spatio-temporal  sidelobes.[8,  9]  The  present  work  develops  a  pressnre  sensitivity 
kernel  that  can  be  used  with  a  variety  of  source  transmissions.  The  sensitivity 
kernels  obtained  for  different  source  transmissions  can  be  used  as  a  performance 
measure  to  study  the  effect  of  medium  sound  speed  finctuations  on  the  received 
signal.  As  an  example,  time-reversal  is  compared  to  a  one-way  transmission,  and 
an  explanation  is  provided  for  the  lower  sensitivity  of  time-reversal.  The  ability  to 
express  the  received  pressnre  perturbation  as  a  function  of  the  sensitivity  kernel 
allows  for  the  development  of  a  regularized  version  of  the  inverse  filter  [10],  with  a 
physically  relevant  regularization.  In  this  paper,  source  transmissions  may  some¬ 
times  be  referred  to  as  source  functions,  or  source  shading  functions  (for  an  array 
of  sonrces)  in  order  to  explicitly  denote  the  time,  freqnency  or  spatial  dependence 
of  the  source  transmission.  This  paper  can  be  seen  as  primarily  providing  a  the¬ 
oretical  framework,  which  is  then  visualized  using  a  simple  acoustic  propagation 
scenario.  Free  space  acoustic  propagation,  followed  by  a  Pekeris  waveguide,  is  used 
to  visualize  various  sensitivity  kernels. 

Wave-theoretic  sensitivity  kernels  were  first  proposed  in  seismology  as  a 
finite-frequency  tomography  alternative  to  the  more  commonly  used  time-of-arrival 
tomography.  [11,  12,  13]  A  key  contribution  was  in  revealing  that  it  is  possible  to 
have  zero  sensitivity  to  sonnd  speed  pertnrbations  along  a  ray  path,  while  having 
non-zero  sensitivities  in  the  region  surronnding  a  ray  path.  More  recently,  travel¬ 
time  sensitivity  kernels  for  ocean-acoustic  tomography  were  developed  and  were 
shown  to  exhibit  similar  properties  to  those  developed  in  seismology.  [14]  These 
kernels  are  cnrrently  being  employed  to  invert  for  sound  speed  perturbations  in  a 
shallow  water  waveguide  using  broadband  aconstic  signals.  [15]  All  of  these  kernels 
make  use  of  a  first-order  Born  approximation,  which  results  in  a  linear  relation 
between  the  sonnd  speed  perturbation  and  either  the  induced  pressure  perturbation 
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in  the  Green’s  function  or  a  travel-time  perturbation.  The  first  Born  approximation 
restricts  the  use  of  these  sensitivity  kernels  to  linear  regimes  where  the  sound 
speed  perturbation  is  small  compared  to  the  background  sound  speed  prohle.  For 
larger  perturbations,  non-linearities  result,  which  can,  to  an  extent,  be  handled  by 
iterative  adjoint  methods.  [16]  As  a  result  of  the  first-order  Born  approximation,  it 
is  assumed  that  the  Green’s  function  used  as  a  source  transmission  in  time-reversal 
represents  an  unperturbed  medium,  and  that  medium  perturbations  are  restricted 
to  a  linear  regime. 

The  contents  of  this  work  are  organized  as  follows.  In  Sec.  II,  pressure 
sensitivity  kernels  are  developed  for  the  general  case  of  an  array  of  hydrophones 
transmitting  to  a  single  receiver  (broadside  transmission).  The  theory  can  be 
easily  extended  to  a  multiple-input  multiple-output  scenario.  Sections  III  and  IV 
show  examples  of  sensitivity  kernels  in  free  space  and  for  a  Pekeris  waveguide, 
respectively.  The  goal  is  to  examine  the  sensitivity  of  time-reversal  versus  that 
of  a  one-way  transmission.  In  a  dispersionless  propagation  environment,  and  in 
the  absence  of  a  transmit  array,  there  is  no  difference  in  sensitivity  between  time- 
reversal  and  a  one-way  transmission.  As  suggested  in  Roux  et  ah, [17]  a  reduction  in 
sensitivity  through  time-reversal  is  only  made  in  the  presence  of  multiple  sources, 
or  multipath  propagation.  Sensitivity  kernels  for  a  source-receiver  path  have  a 
Fresnel  zone-like  interference  pattern,  and  a  reduction  in  sensitivity  for  a  hnite- 
frequency,  hnite-bandwidth  case  is  caused  by  an  averaging  of  several  kernels  that 
correspond  to  multiple  sources  or  multiple  paths.  Time-reversal  results  in  each 
of  the  multiple  paths  arriving  at  the  same  time,  as  a  consequence  of  which  the 
sensitivity  kernels  for  each  of  the  paths  add  coherently  to  give  an  overall  reduction 
in  sensitivity.  This  holds  true  not  just  for  time-reversal,  but  also  for  inverse-hlters, 
or  any  transmission  that  results  in  the  simultaneous  arrival  of  different  paths.  In 
addition  to  the  coherent  cancellation  of  sensitivities,  there  is  also  an  incoherent 
reduction  in  sensitivity  that  comes  from  independent  paths  having  uncorrelated 
sound  speed  perturbations.  The  incoherent  reduction,  however,  is  not  addressed 
in  this  paper. 

In  Sec.  V,  the  sensitivity  kernels  are  examined  from  an  optimization  point 
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of  view,  the  goal  being  to  design  source  functions  that  offer  the  flexibility  of  trad¬ 
ing  mainlobe  width  for  increased  robustness.  Semi-analytic  expressions  for  the 
optimized  signals  are  offered,  which  result  from  a  minimization  of  pressure  per¬ 
turbations  together  with  a  constraint  on  the  focusing  capabilities  of  the  source 
function.  This  approach  leads  to  an  expression  for  the  optimized  source  function 
which  is  an  enhancement  to  the  inverse  £lter.[18]  Sec.  VI  presents  the  results  of  sim¬ 
ulations  in  a  Pekeris  waveguide,  where  the  optimized  source  function  is  compared 
to  time-reversal  and  the  inverse  hlter.  Sensitivity  kernels  for  the  optimized  source 
functions  are  presented,  which  explicitly  use  knowledge  of  fluctuation  statistics  to 
lower  sensitivity  to  medium  perturbations.  Finally,  Sec.  VII  contains  a  discussion 
of  the  results  and  draws  conclusions  from  this  work. 


2.2  Pressure  sensitivity  kernels 

The  pressure  perturbation  to  the  Green’s  function  ^(rlr^)  at  a  location  r 
due  to  a  harmonic  point  source  of  unit  strength  at  location  Vg  can  be  expressed  as 

AG{r\rg)  =  -2u;^  [  G(r|r')G(r>,)^^dr'  (2.1) 

This  is  called  the  first  Born  approximation.  [19]  It  is  also  called  the  single-  or  weak- 
scattering  approximation,  or  wave-held  linearization,  and  is  widely  used  in  connec¬ 
tion  with  scattering  problems.  [20]  It  is  a  linear  relation  between  the  perturbation 
in  the  sound-speed  distribution  and  the  induced  perturbation  in  the  Green’s  func¬ 
tion.  The  integration  volume  in  eq.  2.1  is  all  space  or  the  region  encompassing  the 
source  and  receiver  over  which  sound  speed  perturbations  result  in  perturbations 
to  the  Green’s  function. 

A  key  assumption  upon  which  the  validity  of  the  hrst  order  Born  approxi¬ 
mation  rests  is  that  Ac/c,  integrated  over  the  volume,  be  small.  A  consequence  of 
the  volume  integration  in  eq.  2.1  is  that  the  validity  of  the  hrst  Born  approximation 
in  a  waveguide  will  depend  on  the  path  length.  Longer  path  lengths  will  result 
in  the  linear  approximation  breaking  down  at  smaller  sound  speed  perturbations 
than  for  shorter  path  lengths.  The  ehect  of  non-linearities  was  seen  to  start  to 
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creep  in  for  a  Ac  of  4  ms  ^  with  a  background  sound  speed  of  around  1500  ms 

|21] 

The  Green’s  function  above  is  dehned  in  the  frequency-domain  -  the  cir¬ 
cular  frequency  u  enters  the  Helmholtz  equation  as  a  parameter.  Further,  the 
Green’s  function  depends  on  the  sound-speed  distribution  which  also  enters  the 
Helmholtz  equation.  In  order  to  explicitly  note  these  dependencies,  G  is  written 
in  the  following  as  G{r\rs]uj]c). 

For  a  single  point  source  the  resultant  pressure  held  can  be  written  as 

P(r^;a;;c)  =  G(r^|rs;  cc;  c)S'(rs;  cc)  (2.2) 

where  P{rr',uj-,  c)  is  the  frequency-domain  pressure  at  receiver  location  r^,  from  a 
source  located  at  r^.  G{rr\rs;u;c)  is  the  frequency-domain  Green’s  function  that 
turns  pressure  at  the  source  at  frequency  u  to  pressure  at  the  receiver.  S{rs]uj) 
is  the  signal  at  source  at  frequency  u.  For  an  array  transmitting  to  a  single 
receiver,  the  received  signal  is  given  by  adding  up  all  source  contributions: 

Ns 

P{rr-,u}-,c)  =  ^G{rr\rsi-,  UJ]  c)S{rsi-,  Uj)  (2.3) 

i=l 

where  r**  is  the  location  of  the  source  element,  and  Ng  is  the  number  of  source 
elements. 

In  the  time-domain,  the  received  signal  is  then  given  by 

-^5  poo 

P(r^;t;c)  =  ^/  G{rr\rsi;t  -  t';c)S{rsi;t')dt'  (2.4) 

i=i 

It  is  seen  from  eq.  2.4  that,  in  the  case  of  an  array-to-a-point  transmission,  the  time- 
domain  Green’s  function  has  been  accumulated  over  source  and  time,  weighted  by 
S,  as  opposed  to  a  point-to-point  transmission  where  the  Green’s  function  is  only 
accumulated  over  time.  For  now,  a  heuristic  argument  can  be  made  that  this 
additional  spatial  averaging  could  make  the  received  signal  more  stable.  This  will 
be  shown  more  rigorously  in  Sec.  VI. 

In  a  perturbed  propagation  environment,  the  frequency  domain  perturba¬ 
tion  to  the  received  pressure  can  be  approximated  by  substituting  the  perturbed 
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Green’s  function  into  eq.  2.3 

AP(r^;a;;c)  =  AG(r^|r^^;  a;;  c)5'(r^^; 


N, 


U 


(2.5) 


2=1 


Using  eq.  2.1  for  the  perturbed  Green’s  function 

Ns 


AP(r^;cc;)  =  ^  (^  -  2u;^y  G(r^|r;  a;;  c)G(r|r^i;  o;;  c)^|^)dr  j  P(r*i;  o;)  (2.6) 

Exchanging  the  order  of  integration  and  summation  will  yield  the  effective 
pressure  sensitivity  kernel. 


AP{rr;co)  =  —2uj^ 


C  -^s  A  /  \ 

y  (  S  G{rr\r-,  u;  c)G(r|r^i;  u;  c)  S{r,i;  uj)^dr 


Ns 


'  V  ^  1 

1=1 


'V 


J2Klp{r]U]c)jAc{r)dr 

i=l 

Kap{r-,uj-,  c)Ac{r)dr 


(2.7) 


where  Kp^^v,  u]  c)  is  the  sensitivity  of  the  received  pressure  for  a  receiver  located 
at  Vr  to  a  sound  speed  perturbation  at  location  r,  for  the  single  source  with  index 
i  located  at  r„-. 


Kl  (r;uj;c)  =  -2uj‘^  G{rr\r;u;  c)  G{r\rsi;u;  c)^^^S{rsi;uj) 


c^ir] 


=  Ml(r;uj;c)S{rsi;uj) 


(2.8) 


and  Kap{r;u;  c)  is  the  sensitivity  of  the  received  pressure  to  a  sound  speed  pertur¬ 
bation  at  location  r,  for  an  array  of  sources  transmitting  to  a  single  receiver. 
can  be  recognized  as  part  of  the  integrand  in  eq.  2.1,  i.e.  integrating  Mpp  times 
the  sound  speed  perturbation  over  the  region  between  source  and  receiver  will  give 
the  Born  approximation  to  the  perturbed  Green’s  function. 

Eq.  2.7  is  a  linear  relation  between  the  pressure  perturbation  at  the  receiver 
and  sound  speed  perturbations  in  the  medium,  valid  when  Ac/c  is  sufficiently  small. 

Eqs.  2. 7-2. 8  provide  the  framework  for  studying  the  effect  of  source  trans¬ 
mission  schemes  on  the  received  pressure  perturbation.  Replacing  S{rs]U!)  by 
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c),  will  result  in  an  expression  for  the  pressure  perturbation  for  a  time- 
reversal  scheme.  Replacing  S{rs'iOj)  by  1  /  G{r  j.\r  g-,  uj\  c)  will  result  in  an  expression 
for  the  pressure  perturbation  for  a  scheme  involving  inverse  £lters.[18] 

Finally,  an  additional  way  of  looking  at  eq.  2.7  is  in  terms  of  an  optimization 
problem  where  the  goal  is  to  obtain  a  set  of  source  transmissions  that  will  minimize 
the  pressure  perturbations  at  the  receiver.  This  is  more  fully  explored  in  Sec.  V. 


2.3  Free  space  sensitivity  kernels 

As  a  first  step  towards  visualizing  pressure  sensitivity  kernels,  the  analytical 
framework  developed  in  Sec.  II  is  applied  to  acoustic  propagation  in  free  space. 
The  analytical  framework  for  free  space  closely  follows  that  of  Skarsoulis.  [14] 


2.3.1  Analytical  results 


The  complex  amplitude  of  the  3-D  Green’s  function  between  a  point  source 
located  at  and  receiver  at  in  a  medium  with  uniform  background  sound  speed 
c  has  the  form. 


AkR 


Ggiu)  — 


AttR 


(2.9) 


where  R  =  |r5  —  r^l  is  the  source-receiver  distance,  k  =  u/c. 

Scatterers  in  the  medium  between  source  and  receiver,  each  of  which  con¬ 
tribute  a  sound  speed  perturbation  Ac{x,  y,  z)  will  cause  a  perturbation  to  the 
Green’s  function  (eq.  2.1).  Denoting  the  distance  from  the  source  to  a  scatterer 
(sound  speed  anomaly)  by  Ri  and  from  the  scatterer  to  the  receiver  by  R2,  the 
perturbation  to  the  received  pressure  will  be 


AP{uj)  = 


-‘liKp- 


^ik[R\  (x,y,z)+R2  {x,y,z)'\ 


S{u)Ac{x,y,z)dxdydz  (2.10) 


(47r)2c3y  J  J  Ri{x,y,z)R2{x,y,z) 

From  eq.  2.10  and  eq.  2.8,  the  free  space  sensitivity  kernel  for  a  single 
source-receiver  pair  is 


Kpp{x,y,z;u)  = 


—  luR  Qik[Ri{^^y,z)+R2{x,y,z)] 

{A'kYG’  Ri{x,  y,  z)R2{x,  y,  z] 


Siu) 


(2.11) 
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The  three-dimensional  pressure  sensitivity  kernel  can  be  reduced  to  two-dimensions 
by  integrating  eq.  2.11  over  the  horizontal  cross-range  dimension.  This  assumes 
that  the  sound  speed  anomalies  between  source  and  receiver  do  not  vary  in  the 
horizontal  cross-range  direction.  Next,  the  continuous- valued  kernel  in  eq.  2.11 
can  be  discretized  by  integrating  eq.  2.11  without  the  source  term  over  grid  cells 
(Fig.  2.1),  with  the  sound  speed  anomaly  assumed  to  be  constant  over  each  cell. 
These  two  operations  yield  the  integrated  pressure  sensitivity  for  each  grid  cell 
(eq.  2.8), 


(47r)2c3 


l'Xn+dxl2 
'  Xn  —  dxl2 


^Zm-\-dzl2 
'  Zm-dzl2 


^ik[Ri{x„,y,Zm)+R2ix 

— - — - - r  dzdxdy 

Vi  ^m)R2\^n'j  Vi  ^m) 


(2.12) 


Eq.  2.12  is  numerically  evaluated  in  Section  V. 

Eq.  2.12,  multiplied  by  a  source  function  will  yield  the  point-to-point  sensi¬ 
tivity  kernel  (eq.  2.11),  following  which  the  pressure  perturbation  can  be  calculated 
using  eq.  2.10. 

It  is  helpful  to  look  at  the  time  domain  pressure  sensitivity  equations  for 
the  simple  case  of  free  space  propagation  using  both  a  unity  frequency  domain 
source  function  and  the  phase  conjugated  Green’s  function  as  source  functions. 

For  a  one-way  transmission,  S(u)  =  1  (delta  function  in  time,  centered 
at  t  =  0).  This  will  result  in  a  single  arrival  at  time  t  =  R/c  =  T.  Hence, 
the  time-domain  sensitivity  kernel  for  pressure  at  the  time  of  the  peak  arrival  is 
evaluated  at  time  t  =  R/c.  For  time-reversal,  the  shape  of  the  spectrum  stays  the 
same,  but  the  signal  is  phased  and  shaded  in  such  a  way  that  the  source  function 

—  ikR 

Siu)  ^  .  This  phasing  will  result  in  the  single  arrival  occurring  at  time  f  =  0, 
as  opposed  to  the  one-way  transmission  where  the  arrival  occurs  at  t  =  R/c  =  T. 
Hence,  for  time  reversal,  the  sensitivity  kernel  for  pressure  at  the  time  of  the  peak 
arrival  is  evaluated  at  time  f  =  0.  The  sensitivity  kernel  at  the  time  of  an  arrival 
provides  quantitative  information  about  the  sensitivity  of  the  pressure  amplitude 
of  the  peak  of  the  arrival  to  sound  speed  perturbations  at  any  point  in  the  medium. 
Henceforth,  the  term  time-domain  sensitivity  kernel  is  used  to  denote  time-domain 
pressure  sensitivity  for  a  particular  time,  usually  the  peak  of  an  arrival,  to  sound 
speed  perturbations. 
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Using  the  appropriate  source  function  S{u),  and  applying  the  inverse  Fourier 
transform  to  eq.  2.11  will  give  the  time-domain  sensitivity  kernels  for  point-to-point 
and  time-reversal  transmissions, 


Kpp{i)  — 

Ktr{t)  = 


2 

(47r)2c3i?ii?2 


'  oo 
poo 


(47r)2c3i?ii?2 

1 


(47r)2c3i?ii?2 


5”{Ti  +  T2-t) 

(5"(Ti  +  T2  -  T  -  t) 


(2.13) 


where  S(.)  is  the  Dirac  delta  function  and  '  denotes  its  time  derivative,  Ti  =  Ri/c 
and  T2  =  R2/C.  The  notation  has  been  simplihed  to  drop  the  explicit  dependence 
of  Ri,R2,Ti  and  T2  on  x,y,z.  Eq.  2.13  for  the  time-domain  pressure  sensitiv¬ 
ity  kernels  is  similar  to  what  was  derived  with  respect  to  travel-time  sensitivity 
kernels.  [14] 

Evaluating  the  sensitivity  kernels  at  the  respective  arrival  times  will  result 
in. 

Thus,  in  the  trivial  case  of  a  point-to-point  transmission  in  free  space,  both  time- 
reversal  and  a  one-way  transmission  have  identical  pressure  sensitivities. 

In  a  hnite  frequency,  hnite  bandwidth  case,  the  inhnite  bandwidth  Dirac 
delta  function  is  replaced  by  a  (sine)  form  where  the  argument  of  the  sine 

function  is  given  by  Ti  -|-  T2  —  T.  As  a  result,  the  sensitivity  kernel  in  eq.  2.14 
can  be  approximated  by  the  second  derivative  of  the  sine  function.  This  gives 
rise  to  a  pattern  in  the  range-depth  space  marked  by  lines  over  which  Ti  -|-  T2  is 
constant,  which  is  the  broadband  equivalent  of  the  lines  of  constant  phase  in  a 
single  frequency  Fresnel  interference  pattern. 


2.3.2  Numerical  results 

A  simple  way  to  test  the  sensitivity  kernels  (eq.  2.8)  from  previous  sections 
is  to  conduct  numerical  simulations  involving  acoustic  propagation  in  free  space 


24 


where  there  are  sound  speed  anomalies  6c(x,  z)  in  the  medium  between  source  and 
receiver.  [14] 


Figure  2.1:  Computational  geometry  for  free  space  sensitivity  kernels 

The  scenario  studied  involves  a  two-dimensional  area,  4000  m  in  length  and 
100  m  in  depth.  While  the  theory  developed  here  applies  to  three-dimensional  ge¬ 
ometries,  the  results  are  more  easily  visualized  by  integrating  the  three-dimensional 
pressure  sensitivity  kernel  along  the  horizontal  cross-range  dimension  (eq.  2.12). 
It  was  shown [22]  that  three-dimensional  scattering  effects  become  important  only 
when  the  Fresnel  width  approaches  and  exceeds  the  cross-range  coherence  length 
of  the  sound  speed  perturbations.  This  is  because,  under  such  a  scenario,  out-of- 
plane  scattering  becomes  important,  which  cannot  be  accounted  for  by  2-D  models. 
As  shown  in  the  next  section,  the  Fresnel  width  for  pressure  sensitivity  kernels  in 
the  free  space  scenario  considered  is  approximately  20  m.  The  cross-range  scale 
for  sound  speed  perturbations  was  assumed  to  be  larger  than  30  m,  thus  satisfying 
the  criterion  under  which  a  two-dimensional  example  is  sufficient  to  yield  realistic 
results.  As  a  result,  it  is  possible  to  work  in  two  dimensions  by  using  sensitiv¬ 
ity  kernels  that  are  integrated  in  the  horizontal  cross-range  dimension.  There  are 
10  transmit  transducers  spanning  depths  15  m  to  90  m.  The  receive  transducer 
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is  at  a  range  of  4000  m  and  at  a  depth  of  50  m.  The  background  sound  speed 
is  assumed  to  be  1500  ms“h  The  transducers  operate  at  a  central  frequency  of 
3.5  kHz  (Ac  =  0.4  m)  with  a  bandwidth  of  128  Hz.  These  values  are  chosen  to 
roughly  reflect  those  used  in  recent  ocean  acoustics  experiments.  Fig.  2.1  shows 
the  computational  geometry  employed  to  illustrate  sensitivity  kernels  in  free  space. 
The  size  of  each  grid  cell  is  dx  =  30  m  by  dz  =  2  m.  Sensitivity  kernels  are  nu¬ 
merically  integrated  over  three-dimensions  for  each  grid  cell.  The  size  of  the  grid 
cell  is  representative  of  spatial  coherence  scales  of  sound  speed  perturbations  in 
the  ocean.  [23] 

Broadband  Green’s  functions  were  calculated  from  each  source  to  the  do¬ 
main  and  from  the  receiver  to  the  domain  (a  hnely  sampled  three-dimensional  grid 
cell),  over  the  128  Hz  bandwidth  centered  around  3.5  kHz,  using  a  frequency  step 
of  1  Hz.  The  product  of  the  two  Green’s  functions  is  then  numerically  integrated 
over  each  grid  cell.  This  allows  for  the  calculation  of  the  frequency-domain  sensitiv¬ 
ity  kernel  according  to  eq.  2.12.  An  FFT  is  then  applied  to  the  frequency-domain 
sensitivity  kernels  at  each  grid  point,  which  results  in  a  time-series  of  sensitivity 
for  each  cell.  For  each  cell,  the  value  of  the  sensitivity  time-series  at  the  time  of  an 
unperturbed  arrival  is  the  time-domain  pressure  sensitivity  for  the  arrival.  This 
may  be  expressed  as 

/OO 

M{uj)S{uj)e-^^^duj  (2.15) 

•OO 

For  a  one-way  point-to-point  transmission,  M{u)  =  S{uj)  =  1  and  t  = 

R/c.  For  a  broadside  transmission,  M{u)  =  1  t  =  R/c. 

For  a  time-reversal  mirror,  M{u)  =  J^iLi  ^ppi^)^  t  =  0.  The 

top  left  panel  in  Fig.  2.2  shows  the  time-domain  pressure  sensitivity  kernels  for  the 
main  (only)  arrival,  for  one-way  point-to-point  transmissions  between  a  source  at 
50  m  and  a  receive  element  at  50  m  depth,  at  a  range  of  4  km.  The  structure  of  the 
sensitivity  map  shows  an  approximate  Fresnel-zone  behavior  with  an  oscillatory 
structure.  [14]  The  signihcance  of  the  Fresnel  zones  is  that  within  each  zone,  the 
phase  difference  between  the  interfering  Green’s  functions  is  less  than  tt.  While 
the  Fresnel  zone  is  a  single-frequency  concept,  the  time-domain  analogue,  called 
the  ’zone  of  influence’  [24]  is  a  weighted  average  of  many  single-frequency  Fresnel 


26 


interference  patterns.  As  the  central  frequency  increases,  the  width  of  the  zone 
of  influence  becomes  smaller  as  the  propagation  becomes  more  ray-like.  As  the 
bandwidth  increases,  the  higher  order  Fresnel  zones  of  different  frequencies  begin 
to  destructively  interfere,  and  only  the  hrst  zone  of  influence  is  discernible.  The 
width  of  the  central  lobe  of  the  kernel  is  approximately  the  same  as  that  of  the 
hrst  Fresnel  zone  at  the  central  frequency,  given  by  [25] 

«.(.>. B..1 

where  R  is  the  source-receiver  distance,  /  is  the  central  frequency  of  the  signal, 
and  c  is  the  sound  speed.  Taking  R  =  4000  m,  r  =  2000  m,c  =  1500  m/s  and 
/  =  3500  Hz,  gives  Rp  =  20.7  m,  which  is  consistent  with  the  observed  value  of 
approximately  20  m  in  Fig.  2.2.  As  a;  — oo,  Rri'f')  0,  and  the  wave-theoretic 
sensitivity  kernel  becomes  the  ray-theoretic  kernel.  The  sensitivity  on  the  source- 
receiver  path  is  non- zero,  which  is  different  from  3-D  travel-time  sensitivity  kernels 
which  have  zero  sensitivity  along  the  ray  path.  In  seismic  literature,  [11]  3-D  travel¬ 
time  sensitivity  kernels  are  often  referred  to  as  ’banana-doughnut’  kernels  with  the 
’banana’  shape  arising  from  the  shape  of  the  kernel  as  a  result  of  refraction  of  the 
ray  path,  and  the  ’doughnut’  referring  to  the  region  of  zero  sensitivity  along  the 
ray  path.  The  non-zero  sensitivity  along  the  ray  path  has  also  been  observed  for  2- 
D  travel-time  kernels.  [26]  Because  the  2-D  kernel  is  obtained  from  the  cross-range 
averaging  of  the  3-D  kernel,  the  2-D  kernel  has  non-zero  sensitivity  along  the  ray 
path. 

With  an  understanding  of  the  pressure  sensitivity  kernels  in  place,  it  is  now 
possible  to  study  the  effect  of  time-reversal  on  the  sensitivity  kernel.  Any  reduction 
in  sensitivity  obtained  through  time-reversal  [27]  may  be  thought  of  as  a  two-part 
process.  First,  there  is  the  reduction  in  sensitivity,  if  any,  obtained  through  the 
use  of  a  time-reversed  Green’s  function  as  a  source  function.  Second,  there  is  the 
effect  of  using  an  array  of  sources,  which  involves  a  phased  sum  of  several  one¬ 
way  point-to-point  kernels  such  as  shown  in  Fig.  2.2.  Thus,  when  studying  the 
effect  of  time-reversal  on  sensitivity,  it  is  helpful  to  break  up  the  problem  into  a 
time- reversal  part  followed  by  an  array  gain  part,  or  vice  versa. 
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Figure  2.2:  Time  domain  sensitivity  kernels  (10“®/xPa  s  m~^)  evaluated  at  the 
peak  of  the  arrival  for;  Top  Left  Point-point  one-way,  Top  Right  Array-to-point 
one-way  (broadside)  ,  Bottom  Left  Array-to-point  time-reversal  (time-reversal 
mirror).  Receiver  is  on  the  right  hand  side  of  each  plot,  at  a  depth  of  50  m. 


The  top  right,  and  bottom  panels  in  Fig.  2.2  show  sensitivity,  in  the  time- 
domain,  for  one-way  and  time-reversal  transmissions  involving  an  array  transmit¬ 
ting  to  a  single  receiver.  In  the  case  of  time-reversal,  the  power  of  the  time-reversed 
Green’s  function  has  been  normalized  to  have  the  same  total  power  as  the  uniform 
source  function  used  in  a  one-way  transmission.  Since  the  sensitivity  kernels  for  the 
array-to-point  kernels  are  a  sum  of  individual  point-to-point  kernels,  the  array-to- 
point  kernels  have  been  divided  by  the  number  of  sources.  The  broadside  sensitivity 
kernel  is  computed  by  hrst  summing  the  individual  frequency-domain  sensitivity 
kernels  for  each  source,  followed  by  evaluating  the  broadside  time-domain  sensi¬ 
tivity  kernel  at  the  time  of  arrival  of  a  pulse  from  a  source  at  50  m  depth.  This 


is  because  the  received  signal  for  a  broadside  transmission  consists  of  just  a  single 
arrival  corresponding  to  a  source  at  50  m  depth  since  the  difference  in  range  be¬ 
tween  the  various  sources  and  the  receiver  is  insufficient  to  produce  several  distinct 
direct  arrivals.  If  there  were  indeed  distinct  arrivals,  the  broadside  time  domain 
sensitivity  kernel  would  have  to  be  evaluated  at  each  of  the  arrival  times. 

The  r.m.s  values  of  sensitivities  for  the  various  schemes  are  listed  in  Ta¬ 
ble  2.1,  divided  by  the  r.m.s  value  for  a  one-way  transmission.  The  r.m.s  values 
are  divided  by  the  peak  amplitude  of  the  arrival  for  the  one-way  transmissions,  or 
the  peak  amplitude  of  the  focal  spot  for  time-reversal.  The  statistical  signihcance 
of  using  the  r.m.s  of  the  sensitivity  kernel  follows  from  eq.  2.7.  For  uncorrelated 
perturbations,  the  r.m.s  of  the  sensitivity  map  is  shorthand  for  the  r.m.s  of  the 
expected  pressure  perturbation. 

Point-to-point  time-reversal  is  identical  to  a  one-way  transmission  in  terms 
of  sensitivity  (eq.  2.14).  The  use  of  time-reversal  shows  a  reduction  in  sensitivity 
in  free  space  only  when  used  along  with  an  array.  For  an  array  in  free  space, 
time-reversal  results  in  a  phasing  of  transmissions  from  different  sources  such  that 
the  peak  of  all  arrivals  occurs  at  the  same  time.  The  phased  addition  of  sev¬ 
eral  point-to-point  sensitivity  kernels  results  in  a  net  cancellation  of  the  oscillating 
Fresnel-like  zones  of  influence.  A  given  sound  speed  perturbation  will  be  in  the 
zone  of  influence  of  sensitivity  kernels  for  multiple  paths.  When  these  paths  are 
phased  such  that  they  all  arrive  at  the  same  time,  then  evaluating  the  sensitivity 
at  that  time  corresponds  to  an  average  of  all  the  sensitivity  kernels.  Due  to  over¬ 
lapping  positive  and  negative  regions  of  sensitivity,  there  is  an  overall  reduction 
in  sensitivity.  As  frequencies  get  larger,  there  are  two  effects  that  need  to  be  con¬ 
sidered.  First,  increasing  the  bandwidth  results  in  the  addition  of  kernels  for  each 
frequency,  and  this  smears  out  the  higher  order  zones  of  influence  which  reduces 
the  overall  r.m.s  sensitivity.  At  the  same  time,  there  is  less  overlap  of  sensitivity 
kernels  for  the  multiple  paths.  Second,  increasing  the  central  frequency  results  in 
a  larger  pressure  sensitivity  due  to  larger  phase  shifts.  The  combined  effect  is  that 
time-reversal  involving  higher  frequencies  has  a  greater  sensitivity  to  sound  speed 
perturbations,  as  has  been  experimentally  observed.  [5] 


29 


In  the  example  shown,  the  time-reversal  mirror  reduces  the  sensitivity  of 
the  peak  magnitude  to  about  two-thirds  that  of  a  broadside  transmission.  While 
the  broadside  transmission  in  this  example  also  results  in  a  single  arrival,  the 
alignment  of  the  different  peaks  is  not  perfect,  and  the  sensitivity  of  the  peak  of 
the  arrival  is  larger  than  that  for  time-reversal. 

Table  2.1:  The  r.m.s  values  for  one-way  and  time-reversal  free  space  sensitivity 
kernels,  normalized  to  the  sensitivity  for  a  one-way  transmission 


One-way 

(point-to-point) 

Time-reversal 

(point-to-point) 

Broadside 

(array-to-point) 

Time-reversal  mirror 
(array-to-point) 

1 

1 

0.48 

0.33 

2.4  Pekeris  waveguide  sensitivity  kernels 


The  complex  amplitude  Green’s  function  in  a  two-dimensional  Pekeris  waveg¬ 
uide  for  a  source  located  at  (r*,  Zg)  and  a  receiver  at  (r,  z)  is  given  by  [28] 

M 

G{r,  Z)  =  ^  ^  '4)rn{Zs)'ipm{z)H^Q\krmr)  (2.17) 

^  m=l 


For  mode  m,  krm  is  the  horizontal  wavenumber,  is  the  vertical  mode  shape  and 
Hl^\krmX)  is  the  Hankel  function  of  the  hrst  kind.  Eq.  2.17,  along  with  eq.  2.8, 
will  give  the  frequency-domain  pressure  sensitivity  kernel  for  a  Pekeris  waveguide. 
The  waveguide  geometry  for  which  the  Pekeris  sensitivity  kernels  are  computed  has 
dimensions  identical  to  Fig.2.1,  with  the  surrounding  absorbing  medium  replaced 
by  a  pressure  release  surface  and  an  absorbing  bottom.  The  density  and  sound 
speed  in  the  bottom  are  1800  kg/m^  and  1800  m/s  respectively.  While  it  was 
relatively  straightforward  to  compute  closed  form  sensitivity  kernels  for  free  space 
propagation,  the  same  is  not  true  for  the  Pekeris  waveguide.  However,  it  can  now 
be  expected  that  multipath  propagation  will  result  in  differences  between  time- 
reversal  and  one-way  transmissions,  as  opposed  to  free  space  propagation  where 
the  two  are  identical  for  a  point-to-point  application. 

Frequency-domain  pressure  sensitivity  kernels  are  computed  every  1  Hz 
over  a  bandwidth  of  1  kHz  with  a  central  frequency  of  3.5  kHz.  The  two- 
dimensional  grid  over  which  the  sensitivity  kernels  are  computed  is  sampled  well 
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above  the  spatial  Nyquist  frequency  in  order  to  avoid  spatial  aliasing,  with  a  depth 
step  of  dz  =  0.1  m  and  a  range  step  of  dr  =  10  m.  Unlike  the  3-D  free  space  sen¬ 
sitivity  kernels,  which  were  averaged  over  the  cross-range  dimension,  no  averaging 
is  needed  for  the  2-D  Pekeris  waveguide  example,  since  the  Green’s  function  in 
eq.  2.17  is  specifically  dehned  for  a  2-D  geometry.  Green’s  functions  for  this  exam¬ 
ple  are  on  the  order  of  10“^,  as  a  result  of  which  sensitivity  kernels  will  be  on  the 
order  of  10“^°.  In  order  to  avoid  plotting  these  extremely  small  values.  Green’s 
functions  are  scaled  up  by  10^. 

An  inverse  FFT  of  the  frequency-domain  sensitivity  kernels  in  eq.  2.8  will 
result  in  a  time  series  for  each  grid  point.  This  time  series  is  then  sampled  at 
times  corresponding  to  the  time  of  arrival  of  the  various  multipaths  for  a  pulse 
propagating  in  the  unperturbed  medium.  This  will  result  in  a  number  of  sensitivity 
kernels,  each  of  which  corresponds  to  a  particular  multipath  arrival.  The  resulting 
sensitivity  kernels  are  broadband,  and  centered  around  the  carrier  frequency. 

Fig.  2.3  shows  the  time-domain  sensitivity  kernels  for  a  source  located  at 
a  depth  of  50  m  and  a  receiver  located  at  the  same  depth,  but  1000  m  away 
in  range.  The  same  Fresnel-like  interference  pattern  as  free  space  propagation 
is  seen,  centered  around  each  ray  path.  The  sensitivity  kernel  for  the  bottom 
bounce  is  seen  to  have  the  opposite  polarity  compared  to  the  sensitivity  kernel  for 
the  surface  bounce.  This  is  a  result  of  a  negative  reflection  coefficient  associated 
with  the  surface  bounce.  Another  feature  of  the  sensitivity  kernels  in  Fig.  2.3 
is  the  reduction  in  sensitivity  as  the  launch  angle  gets  steeper.  This  is  a  direct 
consequence  of  Snell’s  law  =  const,  cos  9  =  ^,  ^  ^  0  as  9  ^  n/2),  where 

9  is  the  angle  of  incidence  with  respect  to  the  horizontal.  As  ^  gets  smaller  with 
an  increasing  angle  of  incidence  {9),  a  reduction  in  sensitivity  (^)  follows. 
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Figure  2.3:  Time  domain  sensitivity  kernels  /iPa  s  m~'^)  for  one-way  trans¬ 

mission,  for  a  source  at  50  m  and  receiver  1000  m  away  at  a  depth  of  50  m.  Top: 
Envelope  of  Green’s  function.  A,B,C  are  the  individual  multipaths  corresponding 
to  direct  path,  single  surface/bottom  bounce,  double  bounce  off  the  surface  and 
bottom.  Bottom:  Time-domain  sensitivity  kernel  for  each  arrival  . 
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Figure  2.4:  Time  domain  sensitivity  kernels  jj^Pa  s  m~^)  for  the  peak  of  a 

point-to-point  time-reversal  focal  spot,  for  a  source  at  50  m  and  receiver  1000  m 
away  at  a  depth  of  50  m.  Top:  Pressure  envelope  of  the  focal  spot,  Bottom  Time- 
domain  sensitivity  kernel  for  the  focal  spot,  with  the  time  of  arrival  in  seconds. 


Fig.  2.4  shows  the  time-domain  sensitivity  kernel  for  the  time-reversal  focal 
spot,  obtained  using  the  time-reversed  Green’s  function  in  Fig.  2.3  as  a  source 
transmission.  The  time-reversed  Green’s  function  is  normalized  to  have  unit  power. 
Time-reversal  leads  to  the  different  multipath  arrivals  being  phased  in  a  way  such 
that  they  all  arrive  at  the  same  time,  resulting  in  a  focal  spot.  Similarly,  time- 
reversal  results  in  the  sensitivity  kernel  for  the  focal  spot  being  a  phased  sum  of 
the  sensitivity  kernels  for  the  peak  of  each  arrival.  Table  2.2  lists  r.m.s  sensitivities 
normalized  by  the  r.m.s  sensitivity  of  the  direct  path.  The  peak  amplitude  of  the 
focal  spot  is  different  from  the  peak  amplitude  of  the  multipath  arrivals.  In  order  to 
compare  the  sensitivity  of  the  focal  spot  versus  that  of  the  individual  arrivals,  the 
different  sensitivity  kernels  are  hrst  divided  by  the  peak  pressure  of  the  respective 
arrival. 
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Figure  2.5:  Time  domain  sensitivity  kernels  (10“^*^  /iPa  s  m~‘^)  for  the  peak  of  the 
focal  spot  for  a  time-reversal  mirror.  There  are  10  sources  located  between  30  m 
and  75  m  and  receiver  1000  m  away  at  a  depth  of  50  m.  Top:  Pressure  envelope 
of  focal  spot,  Bottom:  Time-domain  sensitivity  kernel  for  the  focal  spot,  with  the 
time-of-arrival  in  seconds. 


Next,  the  use  of  an  array  (time-reversal  mirror)  is  investigated,  and  the 
resulting  sensitivity  map  is  shown  in  Fig.  2.5.  Ten  sources  are  employed  in  the  same 
conhguration  as  shown  in  Fig.  2.1.  The  resulting  sensitivity  kernels  are  divided 
by  the  number  of  sources,  in  order  to  make  a  comparison  with  Figs.  2. 3-2. 4.  It 
is  seen  that  the  sensitivity  kernel  for  a  time-reversal  mirror  is  a  superposition  of 
sensitivity  kernels  corresponding  to  different  sources,  each  of  which  transmits  a 
time-reversed  Green’s  function. 

As  seen  in  Table  II,  the  use  of  a  time-reversal  mirror  results  in  a  reduction  by 
0.34  in  sensitivity  compared  to  a  single-source  time-reversed  transmission  (point- 
to-point  time-reversal).  This  is  caused  by  the  averaging  of  the  sensitivity  kernels 
corresponding  to  the  different  sources  in  a  manner  similar  to  that  of  the  free  space 
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Table  2.2:  The  r.m.s  sensitivities  for  one-way,  time-reversal  and  a  time-reversal 
mirror  in  a  Pekeris  waveguide,  normalized  to  pressure  sensitivity  on  the  direct 
path. 


One-way 

(point-to-point) 

Time- reversal 
(point-to-point) 

Time- reversal  mirror 
(array-to-point) 

Direct  path 

1 

One  bounce 

0.69 

0.71 

0.37 

Two  bounces 

0.52 

example  in  Sec.  Ill,  in  addition  to  averaging  over  the  multipath  arrivals. 


2.5  Source  function  optimization 

As  mentioned  in  Section  II,  expressing  the  array-to-point  sensitivity  kernel 
with  the  source  transmissions  as  adjustable  parameters  allows  for  optimization  of 
the  source  signal  with  regard  to  the  sensitivity  of  the  received  signal  to  sound  speed 
perturbations.  Other  constraints  on  the  source  signal  can  be  added,  resulting  in 
a  joint  optimization  problem.  In  this  section,  the  optimization  problem  is  posed 
in  terms  of  obtaining  a  spatio-temporal  (depth  and  time)  focal  spot,  while  at 
the  same  time  minimizing  the  effect  of  fluctuations.  The  optimization  is  applied 
in  the  frequency  domain  since  the  pressure  perturbations  caused  by  sound  speed 
fluctuations  (eq.  2.7)  is  more  simply  expressed  in  the  frequency  domain.  As  a 
result  of  the  optimization,  the  quality  of  the  focal  spot  (mainlobe  intensity,  width, 
mainlobe  to  sidelobe  ratio)  may  be  traded  off  for  increased  stability.  [5] 

2.5.1  Minimizing  pressure  perturbations 

In  order  to  state  the  optimization  problem  in  matrix  form,  the  frequency 
domain,  time  domain  and  medium  are  discretized.  To  be  conceptually  simple, 
the  medium  is  parameterized  by  a  two-dimensional  grid  in  range  and  depth.  Thus, 
sound  speed  perturbations  of  the  medium  can  be  specihed  by  the  matrix  Ac,  where 
the  rows  correspond  to  depth  and  columns  to  range.  Similarly,  the  frequency- 
domain  point-to-point  sensitivity  kernel  (eq.  2.7)  over  the  medium  can  be  expressed 
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as  a  matrix  K(a;).  Next,  in  order  to  express  the  integral  in  eq.  2.7  as  the  dot  product 
of  two  vectors,  the  Nz  x  two-dimensional  matrices  Ac  and  K(a;)  are  written 
as  {NzNx)  X  1  vectors  after  stacking  the  columns.  Using  the  vec{-}  notation  to 
denote  a  vectorized  (stacked  columns)  matrix,  eq.  2.7  may  be  written  as 

AP{xr,  Zr^uj)  =  vec{Kap{u)}^vec{Ac}AxAz  (2.18) 

where  the  superscript  (■)^  denotes  matrix  transpose.  The  area  element  (volume 
element  for  3-D)  AxAz  may  be  absorbed  into  Ac  or  K.ap{uj)  so  the  integral  becomes 
a  dot  product.  The  two-dimensional  {Nz  x  N^)  matrix  is  obtained  by  the 
multiplication  of  the  three-dimensional  {Nz  x  N^  x  Ng)  matrix  Mpp  and  the  NgXl 
source  vector  S(a;).  The  multiplication  can  also  be  presented  in  stacked  form, 

vec{Kap{u!)}  =  nec{Mpp(a;)}S(a;)  (2-19) 

where  vec{Wlpp{uj)}  is  the  {NzN^)  x  Ng  matrix  obtained  from  stacking  the  hrst  and 
second  dimension  of  Mpp.  Eq.  2.19  will  yield  the  {NzN^.)  x  1  vector  nec{Kap}  which 
expresses  how  the  arrival  at  the  receiver  is  sensitive  to  medium  perturbations. 

The  received  pressure  perturbation  (eq.  2.18)  can  now  be  expressed  as 

AP{uj)  =  S{ijj)'^  vec{yipp{ijj)}'^  vec{Ac]  (2.20) 

The  portion  of  the  optimization  problem  involving  the  minimization  of  the 
received  pressure  perturbation  can  now  be  conveniently  stated  in  terms  of  a  cost 
function  component 

Ji{uj)  =  <  |APp  > 

=  S(ci;)^  vec{Nipp{uj)}^  <  vec{Ac}vec{Ac}^  >  vec{Nipp{uj)}  S(ci;) 

=  S{u)^  vec{Mpp{uj)}^  Cdc  vec{Mpp{u)}  S(a;) 

=  S(o^)^  C^dc  S{u)  (2.21) 

where  Cdc  is  the  sound  speed  perturbation  covariance  matrix.  The  superscript 
{■)^  denotes  matrix  complex  conjugate  transpose  (Hermitian).  Cmdc  is  the  Ng  x 
Ng  covariance  matrix  of  the  Born-approximation  perturbed  Green’s  functions  at 
frequency  u  between  an  array  of  sources  and  a  single  receiver.  So  element  ij  of 
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Cmdc  is  the  covariance  between  the  pressnre  pertnrbation  at  the  receiver  from  a 
transmission  from  sonrce  i  and  the  pressure  perturbation  at  the  receiver  from  a 
transmission  from  source  j.  Ji  is  not  minimized  separately,  but  is  done  so  as  part 
of  the  joint  optimization  in  order  to  avoid  the  trivial  case  of  S{uj)  =  0. 

The  cost  function  component  (eq.  2.21)  does  not  take  into  account  time  and 
range  shifts  in  the  focal  spot  that  may  occur,  caused  for  example,  by  internal  wave 
induced  sound  speed  fluctuations.  [5]  It  also  requires  each  realization  of  the  focal 
spot  to  be  stable,  without  addressing  the  ensemble  average  of  focal  spots,  which 
was  found  to  have  a  high  stability.  [6] 

The  cost  function  is  computed  independently  for  each  frequency  in  the 
bandwidth  of  interest.  This  has  been  done  in  order  to  be  consistent  with  the 
formulation  of  the  inverse  £lter[18]  which  is  more  conveniently  described  on  a 
single  frequency  basis.  While  the  optimization  can  be  done  across  all  frequencies 
simultaneously,  doing  it  on  a  frequency-by-frequency  basis  is  a  stronger  constraint. 

2.5.2  Obtaining  a  focal  spot 

The  other  part  of  the  optimization  is  the  choice  of  source  function  that  will 
yield  a  focal  spot.  It  was  shown  by  Tauter  et.  al.[18]  that  while  time-reversal  only 
optimizes  the  pressure  amplitude  at  the  focus,  the  inverse  filter  adds  constraints 
on  the  field  around  the  focus,  and  is  thus  able  to  impose  a  desired  sidelobe  level 
or  mainlobe  width.  Since  this  is  a  more  general  optimization  problem,  it  is  clearer 
to  develop  the  optimized  source  function  as  a  robust  version  of  the  inverse  filter 
rather  than  as  a  robust  version  of  time-reversal. 

The  ideal  focal  spot  is  defined  as  a  delta  function  in  depth  and  time  at  the 
range  of  the  receiver, 

u{t,  z)  =  S{z  —  Zr)6{t  —  to)  (2.22) 

In  the  frequency  domain,  u(t,z)  is  written  as  the  Np  x  1  vector  U(ci;),  which  is 
the  desired  focal  spot,  defined  over  Np  depths.  The  problem  now  becomes  one  of 
choosing  S(ci;)  to  produce  a  received  pressure  held  at  a  set  of  Np  depths  that  is  as 
close  to  a  delta  function  as  possible,  given  a  finite  bandwidth.  This  may  be  stated 
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in  the  frequency  domain  as 

UM  =  GMSM  +  AGMSM  +  r  (2.23) 

G  is  a  iVp  X  Ng  matrix  of  unperturbed  Green’s  functions  from  each  source  to  each 
of  the  Np  depths  at  the  range  of  the  receiver.  AG  is  the  (unknown)  perturbation 
to  the  Green’s  function  due  to  the  perturbed  sound  speed  held.  In  the  absence  of 
knowledge  about  AG,  it  is  assumed  that  AGS  is  a  random  variable,  which  can  be 
included  as  part  of  the  error  term,  r.  The  misht  vector  r  is  a  x  1  random  vector 
which  is  assumed  to  be  jointly  normally  (Gaussian)  distributed  with  prior  joint 
probability  density  function  (PDF)  A/'(0,  Cr).  Eq.  2.23  then  reduces  to  a  linear 
model 

U(o;)  =  G(o;)S(o^)  +  r  (2.24) 

While  a  full-depth  representation  is  needed  in  eq.  2.24  in  order  to  obtain 
a  spatio-temporal  focal  spot,  such  a  representation  is  not  needed  in  the  hrst  cost 
function  (eq.  2.21)  where  the  emphasis  is  on  minimizing  received  pressure  sensi¬ 
tivity  to  huctuations  only  at  the  receiver  depth.  Multiple  receiver  depths  can  be 
treated  by  adding  terms  to  the  cost  function  corresponding  to  other  depths. 

A  Bayesian  approach  is  now  adopted  to  hnd  the  optimum  estimator  for 
S(a;)  that  will  satisfy  eq.  2.24.  The  x  1  vector  of  source  excitations  S(a;)  is 
treated  as  an  unknown  random  variable  to  be  estimated.  It  is  assumed  to  be 
jointly  normally  distributed  with  PDF  A/'(0,  Cs)  •  Cs  expresses  constraints  on 
the  size  of  S  and  the  correlation  between  the  signals  from  different  sources.  For 
example,  assuming  Cs  =  /dl  means  that  the  signals  from  different  sources  can  be 
chosen  independently,  but  the  power  of  each  source  is  constrained  to  be  usually  less 
than  2/5.  Assuming  Cr  =  al  means  a  uniformly  distributed  independent  Gaussian 
misht  at  each  depth,  with  variance  a.  This  is  a  convenient  assumption,  but  it  may 
not  be  accurate,  since  the  minimum  width  of  the  focal  spot  is  determined  by  the 
multipaths  available  for  focusing. 

Solving  eq.  2.24  is  equivalent  to  minimizing  the  penalty  function  (where 
the  penalty  is  the  difference  between  the  delta  function  and  the  focal  spot  at  the 
receiver) , 


J2  =  r^Cr'^r  +  S^Cs"^S 


(2.25) 
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J2  may  be  recognized  as  the  negative  of  log  likelihood  function  for  S  given  U.  The 
penalty  function  is  affected  by  the  number  of  depths  at  which  the  misht  vector  is 
enforced,  where  a  large  number  of  closely  spaced  depths  results  in  a  larger  penalty 
function.  Minimizing  J2  with  respect  to  S  will  give  the  maximum  a  posteriori 
estimator,  which  for  the  Gaussian-distributed  case  considered,  is  equivalent  to  the 
minimum  mean  square  estimator,  [29] 

SM  =  (Cs^  +  G^Ck^G)-iG^Ck^U  (2.26) 

The  above  estimator,  with  Cr  =  al  and  Cs  =  /9I  gives  a  generalized  inverse 
hlter[10,  30]  which  does  not  require  G  to  be  full-rank, 

S(o;)  =  (G^G  +  ^I)-^G^U  (2.27) 

P 

If  the  only  concern  is  to  obtain  a  temporal  focal  spot  at  a  single  depth  (in  which 
case  U,  r  and  GG^  are  scalars),  it  is  simpler  to  solve  eq.  2.26  in  its  equivalent 
form  of, 

S(o;)  =  CsG^(GCsG^  +  Cr)'^U  (2.28) 

If  Cs  =  /dl,  eq.  2.28  reduces  to  a  scalar  times  G'^U  (phase  conjugation).  This 
shows  that  for  purely  temporal  focusing,  time-reversal  and  the  inverse  hlter  are 
the  same. 

2.5.3  Joint  optimization 

Adding  the  minimum  sensitivity  cost  function  in  eq.  2.21  to  the  inverse 
hlter  cost  function  in  eq.  2.25  gives 

J3  =  +  S^Cg-'S  +  r^Cr'r  (2.29) 

The  source  function  that  will  minimize  J3  with  respect  to  S  has  a  form  similar  to 
eq.  2.26  except  for  the  added  sensitivity  terms,  and  is  given  by, 

S(o;)  =  (Cs'  +  Cmdc  +  G^Cr-'G)-'G^Cr-'U  (2.30) 

If  Cr  =  al  and  Cg  =  /9I,  then  eq.  2.30  yields, 

S(w)  =  (^I  +  aC„^  +  G"G)-‘G"U 


(2.31) 
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A  few  observations  may  be  made  regarding  the  optimized  source  function  obtained 
in  eq.  2.31.  The  ratio  ^  is  misht  to  transmit  power  ratio.  Setting  a  =  0  (no  error 
tolerated  in  the  focal  spot  criterion)  yields  the  inverse  hlter.  Thus,  eq.  2.31  is  a 
modihcation  of  the  generalized  inverse  hlter,  with  a  stability  term  (aCmdc)  that 
takes  into  account  knowledge  of  the  medium  statistics  and  the  sensitivity  of  the 
Green’s  function  to  sound  speed  perturbations  {Cmdc  =  vec{M}'^  Cdc  nec{M}). 
In  a  medium  with  uncorrelated  perturbations  or  extremely  short  scale  inhomo¬ 
geneities,  pressure  perturbations  will  tend  to  be  uncorrelated  between  sources,  so 
Cmdc  will  tend  to  be  an  identity  matrix,  and  eq.  2.31  tends  to  the  regularized 
inverse  hlter.  Letting  j3  ^  oo  (no  constraint  on  transmit  power)  will  yield  an  al¬ 
ternatively  regularized  inverse  hlter.  Finally,  setting  a  to  large  values  will  result  in 
increased  stability  at  the  cost  of  degradation  of  the  focal  spot  because  a  rehects  the 
allowed  diherence  between  the  observed  focal  spot  and  the  desired  delta-function 
like  focal  spot. 

It  is  useful  to  look  at  the  form  of  eq.  2.31  when  only  a  temporal  focus 
is  desired  (no  focusing  in  depth).  In  such  a  case,  GCr~^G^  is  a  scalar,  and 
an  equivalent  form  of  eq.  2.31,  similar  to  eq.  2.28,  may  be  obtained  by  writing 
=  jl  +  aCmdc-  The  optimized  source  functions  are  then  given  by 

S(o;)  =  PG^(GPG^  +  C,)-^U 

=  APG^U  (2.32) 

It  is  thus  seen  that  the  optimized  solution  for  purely  temporal  focusing  reduces  to 
the  matrix  P  times  the  phase  conjugated  Green’s  function.  If  P  is  a  scaled  identity 
matrix,  then  this  is  again  time-reversal. 

As  seen  in  eq.  2.31,  the  optimized  source  function  requires  rather  complete 
knowledge  of  the  medium  in  order  to  compute  the  sensitivity  kernels,  in  addition 
to  knowledge  of  statistics  of  sound  speed  perturbations.  If  these  were  available, 
then  the  optimized  source  function  can  be  expected  to  perform  better  than  time- 
reversal  but  it  must  be  kept  in  mind  that  an  attractive  feature  of  time-reversal  is 
that  only  the  Green’s  function  between  source  and  receiver  need  be  known. 
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2.6  Optimization  applied  to  a  Pekeris  waveguide 

Studies  of  the  optimized  source  function  with  focusing  in  both  time  and 
depth  are  conducted  for  three  scenarios  involving  different  spatial  structures  of 
sound  speed  perturbations  in  a  Pekeris  waveguide.  For  each  case,  three  trans¬ 
mission  schemes  are  considered:  time-reversal,  the  inverse  hlter  (eq.  2.27),  and 
optimized  source  function  (eq.  2.31).  The  dimensions  of  the  waveguide  and  num¬ 
ber  of  sources  are  the  same  as  in  Sec.  IV.  Sound  speed  perturbations  are  assumed 
to  be  constant  over  a  cell  size  of  10  m  in  range,  1  m  in  depth,  and  a  long  distance  in 
the  horizontal  cross-range  dimension.  The  sound  speed  perturbation  in  each  cell  is 
modeled  as  a  Gaussian-distributed  random  variable.  Changes  in  the  perturbation 
variance  with  depth,  and  horizontal  correlation  scales  of  different  lengths  are  intro¬ 
duced  while  comparing  time-reversal  with  the  optimized  source  shading  function. 
In  the  hrst  case,  there  is  no  correlation  between  perturbations  in  different  cells. 
The  second  case  involves  a  channel  of  zero  variance  in  sound  speed  perturbations 
between  source  and  receiver,  while  the  third  case  has  a  long  decorrelation  scale  in 
both  range  and  depth  for  sound  speed  perturbations. 

For  the  given  medium.  Green’s  functions  were  on  the  order  of  10“^,  resulting 
in  sensitivities  on  the  order  of  10“®.  In  order  to  avoid  plotting  these  extremely  low 
valnes,  the  Green’s  functions  were  scaled  np  by  10"^,  resulting  in  sensitivities  on 
the  order  of  1.  The  scaling  nsed  here  is  different  from  that  in  Sec.  IV  dne  to  the 
fact  that  the  sensitivity  kernels  in  this  section  have  been  snmmed  over  10  m  x  1  m 
cells  in  depth  in  order  to  deal  with  more  manageable  matrices. 

Gost  fnnctions  are  calculated  for  the  various  perturbation  scenarios,  where 
the  cost  fnnction  for  sensitivity,  Ji(a;)  is  compnted  using  eq.  2.21,  and  the  cost  fnnc- 
tion  for  the  sharpness  of  the  focal  spot,  ^2(1^),  is  computed  according  to  eq.  2.25, 
with  the  ideal  focal  spot  specihed  at  11  eqnally  spaced  depths  from  25  m  to  75  m. 
Recognizing  that  the  inverse  hlter  produces  the  best  achievable  focal  spot  for  a 
given  medinm,  the  ideal  focal  spot  (U(c(;)  in  eq.  2.24)  is  chosen  snch  that  the  focal 
spot  cost  fnnction  =  0  for  the  nnregnlarized  inverse  hlter.  The  valne  of  a 

is  now  chosen  snch  that  the  focal  spot  cost  fnnction  for  time-reversal  Jj^(a;)  =  1, 
for  all  freqnencies.  A  valne  of  1  is  chosen  for  the  convenience  of  plotting,  and  can 


41 


be  set  to  a  different  value  in  a  practical  problem.  Next,  (3  can  be  varied  to  yield 
a  focal  spot  cost  function  that  varies  between  0  and  1,  depending  on  the  amount 
of  tradeoff  desired  between  quality  of  focal  spot  and  sensitivity  to  medium  pertur¬ 
bations.  For  all  case  studies  in  this  section,  /?  =  50.  The  values  of  a  and  (3  will 
depend  on  the  specihc  scaling  used  for  the  Green’s  function. 

Sensitivity  kernels  and  source  shading  functions  will  be  the  same  for  time- 
reversal  and  the  inverse  hlter  for  all  three  cases  since  they  are  ignorant  of  the 
presence  or  absence  of  structure  in  the  sound  speed  perturbations.  As  a  result, 
source  shading  functions  and  sensitivity  kernels  for  time-reversal  and  the  inverse 
hlter  are  only  shown  for  the  hrst  case  where  there  is  no  structure  in  sound  speed 
perturbations.  The  received  signals  for  time-reversal  and  the  inverse  hlter  are  also 
the  same  for  diherent  cases  and  hence  are  only  shown  for  the  hrst  case.  Finally,  as 
seen  in  Table  III,  since  the  focal  spot  cost  function  is  very  similar  for  the  inverse 
hlter  and  the  optimized  source  function,  this  is  also  not  shown  for  the  last  two 
case  studies.  Sensitivity  cost  functions  for  the  three  transmission  schemes  (time- 
reversal,  inverse  hlter  and  optimized  source  function),  however,  will  be  diherent  for 
each  case  considered  due  to  diherences  in  the  Cmdc  term  for  each  case  (eq.  2.21). 
The  range  of  values  for  the  sensitivity  cost  function  depends  on  the  specihc  Green’s 
functions  considered  and  the  statistics  of  the  medium.  The  value  of  the  sensitivity 
cost  function  also  depends  on  the  number  of  points  (and  spacing)  at  which  the 
sensitivity  kernel  is  evaluated. 

The  performance  of  time-reversal,  the  inverse  hlter  (eq.  2.27)  and  the  opti¬ 
mized  source  function  under  various  sound  speed  perturbation  scenarios  are  sum¬ 
marized  in  Table  2.3.  The  use  of  an  inverse  hlter  imposes  explicit  constraints  on 
sidelobe  levels  in  both  time  and  depth.  [18]  Since  the  optimized  solution  is  similar 
to  the  inverse  hlter,  in  the  absence  of  noise,  the  focal  spot  obtained  will  have  lower 
sidelobe  levels,  and  hence,  a  smaller  focal  spot  cost  function  (eq.  2.25)  than  time- 
reversal.  The  peak-to-sidelobe  ratios  for  the  inverse  hlter  and  optimized  source 
function  is  about  2  dB  greater  than  that  of  time-reversal.  The  focal  spot  cost 
function  (Table  III)  for  the  inverse  hlter  and  the  optimized  source  function  is  al¬ 
ways  about  90%  less  than  that  of  time-reversal.  However,  although  all  three  source 
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functions  are  normalized  to  have  the  same  source  power,  time-reversal  has  a  higher 
peak  received  power.  This  is  because  time-reversal  is  a  matched  hlter,  which  is  the 
optimal  solution  for  maximizing  received  power,  without  other  constraints. 
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Figure  2.6:  Performance  of  time-reversal,  the  inverse  hlter  and  the  optimized  source 
function.  Sound  speed  perturbations  are  homogeneous,  isotropic  and  uncorrelated 
on  a  10m  x  Im  scale.  Top  left:  Transmit  signals  across  depth.  Top  right: 
Received  signals,  in  time  and  depth,  at  the  focal  spot.  Bottom  left:  Pressure 
sensitivity  kernels  (10“^  /iPa  s  m“^)  for  time-reversal,  the  inverse  hlter  and  the 
optimized  source  function.  Bottom  right:  Cost  functions.  The  focal  spot  cost 
function  for  time-reversal  is  hxed  to  be  a  constant  value  of  1. 


Fig.  2.6  shows  the  results  of  the  numerical  simulation  involving  the  hrst 
scenario  where  the  sound  speed  perturbations  have  the  same  variance  and  are 
uncorrelated  between  grid  points  in  both  dimensions.  The  top  panels  of  Fig.  2.6 
show  the  transmit  and  receive  signals  for  time-reversal  and  the  optimized  source 
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shading  functions.  The  bottom  left  panel  shows  the  sensitivity  kernels  for  each 
scheme  while  the  bottom  right  panel  shows  the  cost  function  metrics  over  a  1  kHz 
bandwidth  for  time-reversal  and  the  optimized  source  function.  The  focal  spot 
cost  function  for  time-reversal  is  hxed  to  be  1,  and  is  not  shown  in  order  to  show 
more  clearly  the  cost  functions  for  the  inverse  hlter  and  the  optimized  source 
function  which  are  on  the  order  of  0.01.  The  lack  of  structure  in  the  sound  speed 
perturbation  results  in  only  a  2.5%  difference  in  sensitivity  between  the  inverse 
hlter  and  time-reversal  (Table  III).  However,  despite  the  lack  of  structure,  the 
sensitivity  cost  function  for  the  optimized  source  function  is  about  35%  less  than 
that  of  time-reversal,  while  the  focal  spot  cost  function  for  both  the  inverse  hlter 
and  optimized  source  function  is  about  90%  less  than  that  of  time-reversal. 


Figure  2.7:  Shape  of  depth-dependent  variance  of  sound  speed  perturbations.  The 
sound  speed  perturbations  obtained  were  used  to  obtain  Fig.  2.8. 

Next,  structure  in  the  sound  speed  perturbations  is  introduced  by  means  of 
a  region  where  the  variance  of  the  sound  speed  perturbations  drops  to  zero.  The 
depth  of  zero  variance  is  placed  at  a  constant  depth  of  52  m.  The  shape  of  the 
depth-dependent  variance  is  shown  in  Fig.  2.7.  The  sound  speed  perturbations, 
however,  are  still  uncorrelated  between  cells.  The  ratio  of  transmit  power  for  the 
source  at  50  m  to  the  average  power  at  other  depths  for  the  optimized  scheme  is 
about  100%  higher  than  that  for  time-reversal,  indicating  a  concentration  of  power 
in  the  unperturbed  region.  For  this  case,  the  inverse  hlter  is  26%  more  sensitive  to 
sound  speed  perturbations  than  time-reversal,  while  the  optimized  source  function 
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Figure  2.8:  Performance  of  time-reversal,  the  inverse  filter  and  the  optimized  source 
function.  Sound  speed  perturbations  are  homogeneous,  isotropic  and  uncorrelated 
on  a  10m  x  Im  scale.  The  variance  of  sound  speed  perturbations  reduces  to  zero 
at  52  m  as  shown  in  Fig.  2.7.  The  source  shading  function  and  sensitivity  kernel 
for  time-reversal  and  the  inverse  filter  are  identical  to  those  in  Fig.  2.6,  and  have 
been  omitted.  Top  left:  Optimized  transmit  signals  across  depth.  Top  right: 
Received  signal  in  time  and  depth,  near  the  focal  spot,  for  the  optimized  source 
function.  Bottom  left:  Optimized  pressure  sensitivity  kernel  (10“"^  /iPa  s 
Bottom  right:  Sensitivity  cost  function.  The  focal  spot  cost  function  is  almost 
identical  to  the  case  in  Fig.  2.6. 


is  about  20%  less  sensitive  compared  to  time- reversal.  The  focal  spot  cost  functions 
are  similar  to  the  previous  case. 
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Figure  2.9:  Performance  of  time-reversal,  the  inverse  filter  and  the  optimized  source 
function.  Sound  speed  perturbations  have  a  long  decorrelation  scale  in  range  (100 
m)  and  depth  (100  m).  The  source  shading  function  and  sensitivity  kernel  for 
time-reversal  and  the  inverse  hlter  are  identical  to  those  in  Fig.  2.6,  and  have  been 
omitted.  Top  left:  Optimized  transmit  signals  across  depth.  Top  right:  Received 
signal,  in  time  and  depth,  near  the  focal  spot,  for  the  optimized  source  function. 
Bottom  left:  Optimized  pressure  sensitivity  kernels  (lO”"^  /iPa  s  m“^).  Bottom 
right:  Sensitivity  cost  function.  The  focal  spot  cost  function  is  almost  identical 
to  the  case  in  Fig.  2.6. 
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The  third  study  (Fig.  2.9)  involves  the  use  of  a  long  (100  m  in  range,  100  m 
in  depth)  spatial  decorrelation  scale  for  the  sound  speed  perturbations.  The  corre¬ 
lation  decays  exponentially  over  100  m  in  range  and  100  m  in  depth.  The  pertur¬ 
bations  are  still  isotropic  and  homogeneous.  This  is  meant  to  represent  low-mode 
internal  waves  in  the  ocean.  The  optimized  source  function  shows  signihcantly 
lower  sensitivity,  with  a  98  %  reduction  (Table  2.3)  compared  to  time-reversal, 
while  the  inverse  hlter  shows  a  52%  reduction  with  respect  to  time-reversal.  At 
the  same  time,  the  quality  of  the  focal  spot  is  similar  to  that  of  the  inverse  hlter 
in  terms  of  peak-to-sidelobe  ratio,  althongh  the  peak  amplitude  has  been  rednced. 
Thus,  the  presence  of  structnre  allows  the  optimized  source  function  to  exploit 
the  knowledge  of  sound  speed  perturbations  to  prodnce  a  focal  spot  that  is  very 
similar  to  that  of  the  inverse  hlter,  bnt  at  the  same  time  considerably  less  sensitive 
to  sonnd  speed  pertnrbations  than  time- reversal  (Table  III). 
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2.7  Conclusions 

A  pressure  sensitivity  kernel  was  developed  for  the  case  where  an  array 
of  sources  transmits  to  a  single  receiver.  This  kernel  is  extendable  to  multiple- 
input  multiple-output  scenarios.  The  Born-approximated  sensitivity  kernel  for  the 
received  pressure  was  used  to  provide  an  explanation  for  why  time-reversal  and  the 
inverse  hlter  are  more  stable  compared  to  a  one-way  transmission.  In  the  absence 
of  multipath,  the  sensitivity  of  a  time-reversed  transmission  for  a  single  source- 
receiver  pair  was  shown  to  be  identical  to  a  one-way  transmission.  A  reduction  in 
pressure  sensitivity  is  only  obtained  through  the  use  of  a  time-reversal  mirror  which 
results  in  the  coherent  addition  of  the  kernels  from  different  sources  as  opposed 
to  an  incoherent  addition  in  the  case  of  a  broadside  transmission.  Time-reversal 
allows  for  the  coherent  addition  of  the  kernels  corresponding  to  different  paths 
arriving  at  the  same  time.  The  presence  of  multipath  has  an  effect  similar  to  a 
time-reversal  mirror,  with  a  lowering  of  sensitivity  resulting  from  the  superposition 
of  the  multipath  kernels. 

The  linear  (first  Born)  approximation  is  only  valid  when  sound  speed  per¬ 
turbations  are  very  small  compared  to  the  background  sound  speed.  Ongoing 
simulations,  beyond  the  scope  of  this  paper,  show  that  the  Born  approximation 
for  a  1  km  by  100  m  waveguide  is  valid  for  sound  speed  perturbations  up  to  about 
0.2  m/s,  for  a  background  sound  speed  of  1500  m/s.  If  only  the  envelope  of  the 
Green’s  function  is  of  interest,  the  range  of  validity  is  extended  until  about  5  m/s. 
This  is  within  the  range  of  most  internal  wave  induced  sound  speed  fluctuations  in 
the  ocean.  Further,  as  a  result  of  the  volume  integral  expression  for  the  perturbed 
Green’s  function,  the  effect  of  longer  path  lengths  is  to  lower  the  range  of  sound 
speed  perturbations  for  which  the  linear  approximation  is  valid. 

The  pressure  sensitivity  kernel  was  used  to  obtain  an  optimized  source  shad¬ 
ing  function  that  can  be  considered  an  environmentally  adaptive  version  of  time- 
reversal  and  the  inverse  hlter.  While  this  optimized  source  function  was  almost 
identical  to  the  inverse  hlter  when  there  is  no  structure  in  the  sound  speed  pertur¬ 
bation  held,  the  gains  are  signihcant  in  the  presence  of  structure,  such  as  a  depth 
of  zero  variance,  or  a  long  decorrelation  length  in  the  sound  speed  perturbations. 
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Chapter  3 


Experimental  Demonstration  of 
Pressure  Sensitivity  Kernels 


Pressure  sensitivity  kernels  were  recently  applied  to  time-reversal  acoustics 
in  an  attempt  to  explain  the  enhanced  stability  of  the  time-reversal  focal  spot 
[Raghukumar  et  al,  J.  Acoust.  Soc.  Am.  124,  98-112  (2008)].  The  theoretical 
framework  developed  was  also  used  to  derive  optimized  source  functions,  closely 
related  to  the  inverse  hlter.  The  use  of  these  optimized  source  functions  results  in 
an  inverse  hlter-like  focal  spot  which  is  more  robust  to  medium  sound  speed  fluc¬ 
tuations  than  both  time-reversal  and  the  inverse  hlter.  The  theory  is  now  applied 
to  experimental  data  gathered  during  the  Focussed  Acoustic  Fields  experiment, 
conducted  in  2005,  north  of  Elba  Island  in  Italy.  Sensitivity  kernels  are  calculated 
using  a  range-independent  sound-speed  prohle,  for  a  geometry  identical  to  that 
used  in  the  experiment,  and  path  sensitivities  are  identihed  to  observed  arrivals. 
The  validity  of  the  kernels  in  tracking  time-evolving  Green’s  functions  is  studied, 
along  with  limitations  as  a  result  of  their  linearized  nature.  An  internal  wave 
model  is  used  to  generate  an  ensemble  of  sound  speed  prohles,  which  are  then  used 
along  with  the  calculated  sensitivity  kernels  to  derive  optimized  source  functions. 
Focal  spots  obtained  using  the  observed  Green’s  functions  with  these  optimized 
source  functions  are  then  compared  to  those  obtained  using  time-reversal  and  the 
inverse-hlter.  It  is  shown  that  these  functions  are  able  to  provide  a  focal  spot 
superior  to  time-reversal  while  being  more  robust  to  sound  speed  fluctuations  than 
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the  inverse  filter  or  time-reversal. 


3.1  Introduction 

The  ability  to  focus  sound  in  a  waveguide  has  been  found  to  be  useful  in 
a  variety  of  applications  such  as  underwater  communications,  target  detection, 
medical  ultrasonics  and  non-destructive  testing  of  materials.  [1,  2,  3]  Time-reversal 
has  been  the  tool  of  choice  in  these  applications  due  its  simplicity  and  robustness 
to  environmental  perturbations.  The  time-reversal  operation  consists  of  first  ac¬ 
quiring  the  Green’s  functions  between  an  array  of  sources  and  a  receiver  at  whose 
location  the  sound  is  to  be  focused.  This  set  of  acquired  Green’s  functions,  when 
time-reversed  and  retransmitted,  produces  a  focal  spot  at  the  receive  location.  Fo¬ 
cusing  is  often  observed  even  after  the  Green’s  functions  have  changed  during  the 
time  that  has  elapsed  between  acquisition  and  retransmission. 

While  the  focal  spot  has  been  shown  to  persist  several  minutes  to  days 
after  the  initial  acquisition  of  the  Green’s  function,  [4,  5]  there  is  a  measurable 
degradation  in  quality  which  eventually  makes  the  re-acquisition  of  the  Green’s 
function  necessary.  In  addition,  the  focal  spot  produced  by  time-reversal  maximizes 
the  energy  at  the  focal  spot  without  constraints  on  sidelobes  in  depth,  range  and 
time.  The  inverse  hlter  has  been  suggested  as  an  alternative  to  time-reversal  which 
can  produce  a  sharper  focal  spot,  [6,  7,  8]  at  the  cost  of  knowledge  required  to  obtain 
additional  Green’s  functions  in  depth,  and  a  reduced  amplitude  of  the  focal  spot. 
However,  the  stability  of  the  focal  spot  with  respect  to  environmental  fluctuations 
is  still  not  constrained.  In  order  to  achieve  enhanced  stability  of  the  focal  spot, 
Kim  et.  al.[5]  proposed  a  modihcation  of  narrowband  time-reversal  which  makes 
the  focus  more  stable,  but  with  considerable  broadening  of  the  focal  spot.  This 
paper  demonstrates  the  practical  use  of  optimized  source  functions  which  are  able 
to  produce  a  sharp  inverse  hlter-like  focal  spot,  while  at  the  same  time  being 
robust  to  environmental  fluctuations,  via  a  sensitivity  kernel-based  regularization 
term  included  in  the  inverse  hlter  formulation. 

Sensitivity  kernels  were  originally  proposed  in  seismology  as  wave-theoretic 
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tomography  tools  as  opposed  to  the  more  commonly  used  ray-theoretic  time-of- 
arrival  tomography.  [9,  10,  11]  In  the  ocean  waveguide,  travel-time  sensitivity 
kernels [12]  were  shown  to  have  the  same  properties  as  those  developed  in  seismol¬ 
ogy,  and  were  used  to  invert  for  sound  speed  perturbations  in  deep  water.  Pressure 
sensitivity  kernels  are  currently  being  employed  as  a  means  of  doing  sound  speed 
inversions  in  shallow  water.  [13] 

Recently,  pressure  sensitivity  kernels  were  used  a  means  of  explaining  the 
remarkable  stability  of  the  time- reversal  focus.  [14]  It  was  shown  that  the  time- 
reversal  operation,  which  causes  the  various  multipath  arrivals  in  a  waveguide  to 
arrive  at  the  same  time,  also  results  in  the  pressure  sensitivity  kernel  for  each  ar¬ 
rival  to  overlap  with  the  kernels  for  other  arrivals,  resulting  in  a  net  reduction  in 
sensitivity  for  time-reversal  compared  to  a  one-way  transmission.  An  offshoot  of 
this  work  was  the  formulation  of  optimized  source  functions  using  pressure  sensitiv¬ 
ity  kernels  that  were  based  on  the  inverse  hlter.  This  paper  demonstrates  the  use 
of  these  optimized  source  functions  using  experimental  acoustic  and  environmen¬ 
tal  data  gathered  during  the  Focussed  Acoustic  Fields  2005  (FAF05)  experiment, 
conducted  off  of  Elba  Island,  Italy  in  July  2005.  The  focal  spot  obtained  using 
these  source  functions  is  seen  to  have  a  peak-to-sidelobe  level  in  between  that  of 
time-reversal  and  the  inverse-filter,  but  a  greater  stability  to  environmental  per¬ 
turbations  than  both  schemes.  The  use  of  optimized  source  functions  requires 
knowledge  of  the  covariance  matrix  for  sound  speed  perturbations  and  sensitiv¬ 
ity  kernels.  This  scheme,  therefore,  should  be  viewed  as  being  synergistic  with  a 
shallow-water  observation  program,  which  requires  similar  knowledge. 

This  paper  is  organized  as  follows.  Sec.  II  briefly  describes  the  Focused 
Acoustic  Fields  experiment.  More  comprehensive  details  may  be  found  in  the 
paper  by  Song  et.  ah  [15]  Acoustic  data  from  this  experiment  is  used  to  compare 
the  time-reversal  focal  spot  with  that  obtained  by  the  optimized  source  functions. 
The  propagation  model  used  to  compute  sensitivity  kernels  is  also  validated  against 
acoustic  data.  Environmental  data  gathered  during  the  experiment  is  used  to 
validate  the  internal  wave  model,  which  in  turn  is  used  to  generate  range-dependent 
sound  speed  prohles.  Sec.  Ill  presents  a  review  of  pressure  sensitivity  kernels 
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previously  derived [14],  and  explores  the  limits  of  the  linear  assumption  used  to 
derive  the  kernels.  The  ability  of  the  sensitivity  kernels  to  track  a  time-evolving 
Green’s  function  is  shown  for  time-scales  on  the  order  of  a  few  minutes,  after 
which  the  linear  assumption  no  longer  holds.  Sec.  IV  discusses  the  optimized 
source  functions  and  compares  the  focal  spot  properties  with  experimental  time- 
reversal  data.  It  is  shown  that  focal  spot  quality  can  be  traded  off  for  robustness 
to  environmental  fluctuations,  and  that  the  optimized  source  functions  can  obtain 
an  inverse  hlter-like  focal  spot,  while  being  robust  to  sound  speed  perturbations. 
Finally,  Sec.  V  concludes  the  chapter  followed  by  an  appendix  which  provides 
details  on  the  internal  wave  model  used. 


3.2  Experiment  overview 

The  FAF05  experiment  was  conducted  north  of  Elba  Island  off  the  west 
coast  of  Italy  in  July  2005. [16,  4,  15]  The  experiment  layout  is  shown  in  Fig.  3.1. 
Acoustic  transmissions  were  made  using  a  29-element  Source  Receive  Array  (SRA), 
and  received  over  a  32-element  Vertical  Receive  Array  (VRA).  The  SRA  spanned 
water  depths  34  m  to  112  m,  with  an  array  element  spacing  of  2.786  m,  while  the 
VRA  spanned  water  depths  48  m  to  110  m,  with  an  array  element  spacing  of  2  m. 
The  water  depth  in  the  vicinity  of  the  experiment  was  121  m,  and  the  arrays  were 
4.071  km  apart  in  range.  As  seen  in  Fig.  3.1,  various  environmental  measurements 
were  made  during  the  course  of  the  experiment,  such  as  surface  waves  using  a 
waverider  buoy  (WR),  temperature  prohles  using  two  thermistor  chains  (TC’s), 
and  current  measurements  using  an  acoustic  doppler  current  meter  (ADCP).  The 
thermistor  chains  had  11  elements  each  and  spanned  water  depths  18.5  m  to  76.5 
m,  with  a  sampling  period  of  10  minutes.  In  addition,  conductivity-temperature- 
depth  (CTD)  measurements  were  made  approximately  every  30  minutes.  In  this 
paper,  only  the  CTD  and  TC  data  are  analyzed. 

Round-robin  pseudo-pulse  transmissions  were  made  every  20  seconds  over  a 
period  of  9  hours  on  2005,  Day  197.  These  transmissions,  meant  to  measure  Green’s 
functions  between  the  source  and  receive  arrays,  consisted  of  1-second  chirps  with 
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a  bandwidth  of  1.4  kHz,  and  a  center  frequency  of  3.5  kHz.  The  received  signals 
on  the  VRA  were  then  match-hltered  to  yield  Green’s  function  estimates. 


Figure  3.1:  Locations  of  various  moorings  during  the  FAF05  experiment.  In  green 
at  the  bottom  is  Elba  Island. 

Fig.  3.2  shows  the  temperature,  salinity  and  sound  speed  prohles  obtained 
from  a  CTD  cast  that  was  made  at  the  beginning  of  the  transmissions,  close  to 
the  VRA.  The  sound  speed  prohle  is  downward  refracting,  with  the  thermocline 
located  around  30  m.  Temperature  measurements  made  by  the  thermistor  chains 
at  the  same  time  are  also  plotted.  While  there  is  a  close  agreement  with  between 
the  CTD  and  the  nearby  TC2,  the  TCI  measurement  is  a  little  different,  show- 
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ing  range  dependence.  Fig.  3.3  shows  the  temperature  fluctuations  measured  on 
thermistor  chain  TC2  during  the  FAF  experiment.  Temperature  fluctuations  are 
plotted  as  Altered  above  and  below  the  inertial  frequency.  Fluctuations  above  the 
inertial  frequency  are  shown  over  a  period  of  24  hours,  while  sub-inertial  frequency 
fluctuations  are  shown  over  12  days.  It  appears  that  there  is  some  internal  wave 
activity  between  the  30-40  m  depth,  with  decreasing  activity  towards  the  surface 
and  bottom.  Also  included  are  the  r.m.s  plots  for  the  high-  and  low-pass  Altered 
temperature  chain  data.  The  shorter  time  scale  of  the  high-pass  Altered  data  al¬ 
lows  for  a  closer  look  at  the  internal  wave  energy  structure  around  the  thermocline. 
The  r.m.s  of  the  temperature  data  is  used  to  constrain  the  energy  parameter  of 
the  internal  wave  model  in  Sec.  III. 


Figure  3.2:  Temperature,  salinity  and  sound  speed  profile.  Crosses  on  tempera¬ 
ture  plot  denote  measurements  made  by  thermistor  chain  TC2,  and  circles  denote 
measurements  made  by  thermistor  chain  TCI,  at  the  same  time.  Crosses  and  dots 
on  the  sound  speed  plot  denote  and  SRA  and  VRA  hydrophone  depths. 


Fig.  3.4  shows  the  evolution,  over  10  minutes,  of  two  observed  Green’s 
function  estimates.  The  first  is  between  a  source  at  112.6  m  and  receiver  at  110  m, 
located  close  to  the  bottom,  while  the  second  is  between  a  source  at  112.6  m,  near 
the  bottom  and  a  receiver  at  48  m,  near  the  surface.  While  both  Green’s  functions 
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Figure  3.3:  Temperature  anomalies  (with  respect  to  mean)  derived  from  ther¬ 
mistor  chains  TCI  (blue  solid)and  TC2  (green  dashed).  The  upper  panel  shows 
fluctuations  above  the  inertial  frequency  over  1  day  while  the  lower  panel  shows 
sub-inertial  frequency  temperature  fluctuations  over  two  weeks.  Also  included  are 
the  r.m.s  plots  for  the  high-  and  low-pass  filtered  temperature  chain  data. 


appear  to  be  relatively  stable,  the  right  panel  shows  more  dynamic  variation  over 
the  10  minute  interval,  due  to  its  location  in  a  region  of  internal-wave  induced 
sound  speed  fluctuations. 


3.3  Environmental  modeling 

In  order  to  compute  the  optimized  source  functions,  it  is  necessary  to  cal¬ 
culate  pressure  sensitivity  kernels  and  estimate  the  covariance  matrix  for  sound 
speed  perturbations.  Using  the  Born  approximations,  sensitivity  kernels  can  be 
derived  from  the  Green’s  functions  between  the  source  and  points  in  the  medium, 
and  between  the  same  points  and  the  receiver.  [14]  Experimental  data  typically  only 
have  available  the  source-receiver  Green’s  functions.  Green’s  functions  are  hence 
computed  via  computer  simulations  and  tested  by  comparison  with  those  that  were 
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Receiver  depth:  110  m 
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Figure  3.4:  Temporal  evolution  of  FAF05  Green’s  functions,  for  a  source  at  112.6 

m. 

experimentally  available.  A  parabolic  equation  model  is  used  to  simulate  acoustic 
propagation  in  a  range- dependent  environment,  and  a  Garrett-Munk-based  inter¬ 
nal  wave  model  is  used  to  generate  time-varying  sections  of  sound  speed  profiles. 

[5] 

To  estimate  the  kernels,  a  single  sound  speed  profile  from  a  GTD  cast 
(Fig.  3.2)  taken  at  the  beginning  of  the  transmissions  was  used  in  a  split-step 
Fade  parabolic  equation  model  (RAM)  [17]  was  used  to  model  sound  propagation 
through  a  range-independent  environment.  [18]  The  parabolic  equation  had  a 
delta-function  starter,  range  and  depth  steps  of  1  m  and  0.05  m  respectively,  and 
6  Fade  coefficients,  which  gives  negligible  phase  errors  over  a  75°  angle  range. 
Green’s  functions  were  computed  every  8  Hz  over  the  2800  -  4200  Hz  bandwidth. 
Fig.  3.5  shows  the  waveguide  parameters  employed  for  the  simulation  of  acoustic 
propagation  for  the  FAF05  environment. 

The  data-derived  Green’s  function  is  obtained  by  match-filtering  the  recep¬ 
tions  using  a  2800  -  4200  Hz  Hanning-windowed  chirp  with  a  roll-off  factor  of  0.2. 
In  order  to  make  a  fair  comparison  with  the  data,  the  same  match-hltering  method 
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Figure  3.5:  Waveguide  parameters  used  in  the  acoustic  propagation  model. 

was  applied  to  the  model-derived  Green’s  functions. 

Fig.  3.6  shows  the  full-depth  Green’s  function  on  the  VRA  obtained  using  a 
source  at  112.6  m.  The  same  source  depth  is  used  for  a  number  of  plots  in  this  pa¬ 
per.  This  is  the  deepest  source  depth  available  in  the  data,  and  gives  the  maximum 
stability  when  transmitting  to  a  deep  receiver.  The  upper  panel  shows  the  experi¬ 
mentally  obtained  Green’s  function,  while  the  lower  two  panels  show  that  obtained 
using  the  parabolic  equation,  for  a  range  independent  and  range  dependent  sound 
speed  prohle.  Range-dependent  sound  speed  prohles  were  obtained  using  the  in¬ 
ternal  wave  model,  described  in  the  following  paragraphs.  The  range-dependent 
calculation  does  not  differ  greatly  from  the  range-independent  calculation  because 
the  horizontal  decorrelation  scale  is  about  4  km,  not  much  smaller  than  the  range 
itself.  Almost  all  the  paths  in  the  experimental  Green’s  function  are  resolved  by 
the  parabolic  equation.  The  close  visual  agreement  between  the  two  allowed  for 
the  use  of  an  experimentally  obtained  range- independent  sound  speed  prohle  in 
order  to  compute  the  sensitivity  kernels. 

Internal  waves  are  a  major  source  of  sound  speed  huctuations,  particularly 
in  the  absence  of  any  signihcant  surface  wave  activity,  as  was  the  case  during 
the  course  of  the  FAF05  experiment.  In  addition,  the  downward-refracting  sound 
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Figure  3.6:  Green’s  functions,  model  v/s  data,  for  a  source  at  112.6  m 

speed  profiles  encountered  (Fig.  3.2)  minimize  the  effects  of  surface  waves  on  the 
energetic  arrivals.  A  modihed  Garrett-Munk  internal  wave  model,  [5,  19]  was  used 
to  calculate  vertical  displacements,  from  which  sound  speed  fluctuations  were  cal¬ 
culated  according  to, 

5c{r,  z,  t)  =  Cp[z  -  r]{r,  z,  t)]  -  Cp{z)  (3.1) 

where  Cp{z)  is  the  potential  sound  speed  and  ri{r,z,t)  is  the  vertical  displacement 
of  internal  waves.  The  method  used  to  calculate  the  vertical  displacements  due  to 
the  internal  waves,  identical  to  that  followed  by  Kim  et.  al.  [5],  is  summarized  in  the 
Appendix.  Model  parameters  such  as  the  internal  wave  spectrum,  vertical  mode 
shapes,  mode  amplitudes  and  variance  of  temperature  fluctuations  are  validated 
against  the  thermistor  chain  data. 

Fig.  3.7  shows  a  representative  buoyancy  prohle  calculated  from  a  GTD  cast 
whose  details  are  shown  in  Fig.  3.2.  The  internal  wave  spectrum  calculated  using 
this  buoyancy  profile  is  also  shown,  along  with  the  power  spectrum  of  the  temper¬ 
ature  time  series.  The  energy  parameter  for  the  internal  wave  model  (Appendix) 
is  constrained  by  choosing  a  value  such  that  the  peak  standard  deviation  of  model 
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Figure  3.7;  Top  left:  Buoyancy  profile  calculated  from  CTD  data,  Top  right: 
G-M  model  spectrum  (dashed)  v/s  thermistor  chain  spectrum  (solid),  Bottom 
left:  Standard  deviation  of  model-based  and  thermistor  chain  temperature  fluctu¬ 
ations,  Bottom  right:  G-M  model  mode  amplitudes  v/s  temperature  EOF  mode 
amplitudes 


fluctuations  matches  that  of  the  temperature  fluctuations,  high-pass  hltered  above 
the  inertial  frequency.  Even  so,  the  energy  parameter  that  is  used  is  probably  too 
low  due  to  the  sharp  thermocline  in  the  single  prohle.  Undoing  the  effect  of  these 
sharp  transitions  is  currently  work  in  progress.  Mode  amplitudes  as  used  by  the 
model  are  compared  to  temperature  chain  data  by  projecting  the  temperature  data 
onto  its  EOF  modes  to  yield  the  mode  amplitudes  for  the  temperature  time  series. 
Mode  shapes  are  also  compared  to  the  temperature  data,  and  the  hrst  four  mode 
shapes  are  shown  in  Fig.  3.8.  The  hrst  mode  two  mode  shapes,  which  contribute 
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about  75%  (Fig.  3.7)  of  the  total  internal  wave  energy,  show  a  match  between  the 
model  and  data.  Details  of  the  method  used  to  calculate  mode  shapes  and  mode 
amplitudes  may  be  found  in  the  Appendix. 


Amplitude 


Figure  3.8:  Model  mode  shapes  (dashed)  v/s  data-derived  EOF  mode  shapes 
(solid) 


The  hnal  model  test  involves  combining  the  internal  wave  model  and  acous¬ 
tic  propagation  code,  and  comparing  the  behavior  of  the  time-reversal  focal  spot 
between  theory  and  observation.  Ten  source  transducers  were  chosen,  spaced  5.6 
m  apart,  spanning  the  same  depths  as  ten  similarly-spaced  SRA  transducers.  The 
internal  wave  model  was  used  to  generate  2  days  of  range-dependent  sound  speed 
prohles.  The  parabolic  equation  model  was  then  implemented  for  the  hrst  10  min¬ 
utes  of  sound  speed  prohles.  The  hrst  set  of  Green’s  functions  between  the  source 
array  and  the  receiver  is  time-reversed  and  used  as  the  source  shading  function  in 
order  to  compute  the  time-evolving  focal  spot  for  the  model.  Next,  focal  spots  for 
the  FAF05  data  are  computed  by  taking  10  minutes  of  the  Green’s  function  data, 
time- reversing  the  hrst  set,  and  match-hltering  with  the  rest  to  give  time-evolving 
focal  spots.  Fig.  3.9  shows  the  evolution  of  the  peak  pressure  for  the  time-reversal 
focal  spot  at  four  diherent  depths  in  the  water  column.  The  vertical  bars  represent 


65 


90%  confidence  intervals  obtained  using  9  different  lO-minute  segments  of  Green’s 
functions.  While  the  model  is  able  to  capture  the  initial  slope  in  the  evolution  of 
the  focal  spot,  some  observed  changes  in  slope  are  not  reproduced.  This  limitation 
is  particularly  relevant  in  the  context  of  pressure  sensitivity  kernels,  which  will  be 
discussed  in  the  next  section.  The  ability  to  model  the  slope,  however,  is  a  test  of 
the  validity  of  the  energy  parameter  used  in  the  internal  wave  model.  For  example, 
a  larger  energy  parameter  should  increase  the  rate  at  which  the  focal  spot  pressure 
varies.  This  method  can  also  be  used  as  a  means  to  constrain  the  energy  parameter 
of  an  internal  wave  model. 
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Figure  3.9:  A  comparison  of  the  evolution  of  the  time-reversal  focal  spot  peak 
amplitude  for  the  model  and  data. 


3.4  Pressure  sensitivity  kernels 

3.4.1  Background 

This  section  essentially  summarizes  the  derivation  of  pressure  sensitivity 
kernels  as  presented  in  Skarsoulis  et.  al.[12] 
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The  pressure  perturbation  to  the  Green’s  function  G(r|rg)  between  at  a 
location  r  due  to  a  harmonic  point  source  of  unit  strength  located  at  Vg  can  be 
expressed  as 

AG(r|r,)  =  -20;^  [  G(r|r')G(r>,)^^dr'  (3.2) 

where  u  the  circular  frequency  of  the  source,  and  c(r)  is  the  sound-speed  distribu¬ 
tion. 

This  is  called  the  first  Born  approximation.  [20]  It  is  also  called  the  single- 
or  weak-scattering  approximation,  or  wave-held  linearization,  and  is  widely  used 
in  connection  with  scattering  problems.  [21]  It  is  a  linear  relation  between  the  per¬ 
turbation  in  the  sound-speed  distribution  and  the  induced  perturbation  in  the 
Green’s  function.  The  integration  volume  in  eq.  3.2  is  all  space  or  the  region  en¬ 
compassing  the  source  and  receiver  over  which  sound  speed  perturbations  result 
in  perturbations  to  the  Green’s  function. 

A  key  assumption  upon  which  the  validity  of  the  hrst  order  Born  approxi¬ 
mation  rests  is  that  Ac/c,  integrated  over  the  volume,  be  small.  The  limits  of  the 
hrst  Born  approximation  in  connection  with  experimental  data  are  explored  in  the 
next  sub-section. 

For  a  single  point  source  the  resultant  pressure  held  can  be  written  as 

P{rr;u;c)  =  G{rr\rs;u;c)S{rs;u)  (3.3) 

where  P{rr',uj-,  c)  is  the  frequency-domain  pressure  at  receiver  location  r^,  from  a 
source  located  at  Vg.  G{rr\rg;u;c)  is  the  frequency-domain  Green’s  function  that 
turns  pressure  at  the  source  at  frequency  u  to  pressure  at  the  receiver.  S{rg]u) 
is  the  signal  at  source  Vg  at  frequency  u.  For  an  array  transmitting  to  a  single 
receiver,  the  received  signal  is  given  by  adding  up  all  source  contributions: 

Ns 

P{rr;uj;c)  =  '^G{rr\rgi;uj;c)S{rgi;uj)  (3.4) 

i=l 

where  Vgi  is  the  location  of  the  source  element,  and  Ng  is  the  number  of  source 
elements. 

In  a  perturbed  propagation  environment,  the  frequency  domain  perturba¬ 
tion  to  the  received  pressure  can  be  approximated  by  substituting  the  perturbed 
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Green’s  function  into  eq.  3.4 

Ns 

AP(r^;a;;c)  =  ^  AG(r^|r^j;  a;;  c)5'(r^j;  c^;)  (3.5) 

i=l 

Using  eq.  3.2  for  the  perturbed  Green’s  function  yields  and  expression  for  the 
pressure  perturbation  at  the  receiver  for  an  arbitrary  source  function, 

f  Ac(r) 

AP{rr-,ui)  =  /  }  -2uj‘^  G{rr\r-,ui-,c)  G{r\rsi-,ui-,c)^^^S{rsi]uj)dr 

Jv^  c3(r) 

=  f  Kap{r;u;c)Ac{r)dr  (3.6) 

Jv 

where  K{r]U]c)  is  the  sensitivity  of  the  received  pressure  to  a  sound  speed  per¬ 
turbation  at  location  r,  for  an  array  of  sources  transmitting  to  a  single  receiver. 

Eq.  3.6  is  a  linear  relation  between  the  pressure  perturbation  at  the  receiver 
and  sound  speed  perturbations  in  the  medium,  valid  when  Ac/c  is  sufficiently  small. 
It  provides  the  framework  for  studying  the  effect  of  a  source  transmission  scheme 
on  the  received  pressure  perturbation.  Replacing  S{rs]uj)  by  G*{rr\rs]uj]c),  will 
result  in  an  expression  for  the  pressure  perturbation  for  time-reversal.  Replacing 
S{rs;uj)  by  l/G{rr\rs',uj-,c)  will  result  in  an  expression  for  the  pressure  pertur¬ 
bation  for  the  inverse  £lter.[22]  In  Sec.  V,  eq.  3.6  is  examined  in  terms  of  an 
optimization  problem  where  the  goal  is  to  obtain  a  set  of  source  transmissions 
that  will  minimize  pressure  perturbations  at  the  receiver. 

3.4.2  Experimental  demonstration 

In  this  section,  the  validity  of  pressure  sensitivity  kernels,  briefly  summa¬ 
rized  in  the  previous  section  and  the  references  contained  therein,  is  tested  using 
experimental  data  collected  during  the  FAF05  experiment. 

The  same  parabolic  equation  code  mentioned  in  Sec.  III-A  is  used  to  cal¬ 
culate  the  pressure  sensitivity  kernel  (eq.  3.6)  for  a  given  source-receiver  pair.  For 
every  source-receiver  pair,  two  forward  runs  of  the  code  are  required.  The  first  with 
the  source  placed  at  the  source  depth,  propagated  out  to  the  range  of  the  receiver, 
and  the  second  with  the  source  placed  at  the  depth  of  the  receiver,  propagated  out 
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to  the  same  range.  The  results  of  the  second  run,  flipped  in  the  range  dimension 
is  multiplied  at  each  range  and  depth  step  with  the  result  of  the  hrst  run  to  yield 
the  pressure  sensitivity  kernel  at  each  range  and  depth  step,  for  each  frequency 
(eq.  3.6).  The  background  sound  speed  prohle  at  which  the  sensitivity  kernels  are 
calculated  is  obtained  from  one  of  the  CTD  casts,  taken  at  the  beginning  of  the 
transmissions  (Fig.  3.2).  Sensitivity  kernels  are  calculated  for  a  range-independent 
sound  speed  prohle,  on  top  of  which  range-dependent  sound  speed  perturbations 
are  imposed  in  order  to  compute  Green’s  function  perturbations. 

Fig.  3.10  shows  the  pressure  sensitivity  kernels,  in  the  time-domain,  for  a 
source  at  a  depth  of  112.6  m  and  a  receiver  at  a  depth  of  110  m.  These  depths 
correspond  to  the  depths  of  the  deepest  element  of  the  SRA  and  VRA  respectively. 
The  Green’s  function  is  shown  on  the  upper  panel,  with  various  paths  identihed  as 
A-F.  FAF05  observations  of  the  Green’s  functions  for  the  same  source  and  receiver 
depths  can  be  seen  in  Fig.  3.4.  Time-domain  sensitivity  kernels  are  calculated  by 
taking  the  inverse  Fourier  transform  of  the  frequency  domain  sensitivity  kernels 
and  then  evaluating  the  resulting  time-series  at  the  time-of-arrival  of  the  various 
paths. 

The  pressure  sensitivity  kernels  in  Fig.  3.10  show  an  oscillatory  structure 
around  the  ray  paths  which  is  associated  with  Fresnel-zone  behavior [12],  aris¬ 
ing  from  the  interference  of  two  plane  waves  (eq.  3.6). [23]  While  the  Fresnel- 
zone  is  a  single  frequency  concept,  the  broadband  analogue,  called  the  ”zone- 
of-influence”  [24] ,  is  a  superposition  of  several  single-frequency  Fresnel  zones.  The 
effect  of  bandwidth  then  is  a  decaying  of  the  oscillatory  structure  away  from  the 
hrst  zone,  with  increased  bandwidth  resulting  in  a  faster  decay.  The  effect  of  in¬ 
creasing  central  frequency  is  a  narrowing  of  the  width  of  each  zone-of-inhuence. 
The  combined  effect  of  increasing  bandwidth  and  central  frequency  is  that  the 
wave-theoretic  sensitivity  kernel  starts  to  approach  the  ray-theoretic  kernel. 

Next,  the  limits  of  the  linear  hrst  Born  approximation  used  to  derive  the 
pressure  sensitivity  kernels  (eq.  3.2),  are  tested  numerically  for  the  sound  speed 
prohles  observed  in  the  experiment.  The  linearity  test  involves  computing  a  per¬ 
turbed  Green’s  function  via  a  forward  model  run,  and  then  comparing  with  the 
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Figure  3.10:  Sensitivity  kernels  for  FAF05  paths,  for  a  source  at  112.6  m,  and 
receiver  at  110  m.  The  topmost  panel  shows  the  specihc  Green’s  function  for 
whose  paths  sensitivity  kernels  are  shown  in  the  lower  panels.  Units  of  sensitivity 
are  fiPa  s  m~^.  FAF05  Green’s  function  observations  can  also  be  seen  in  Fig.  3.4. 

estimate  obtained  using  the  sensitivity  kernel.  The  estimate  produced  using  the 
sensitivity  kernel  is  given  by 

^kernel  (^)  ^orig  (i)  +  X]  K{r,  z,t)6c{r,  z)  (3.7) 

r,z 

where  Gorigif)  is  the  unperturbed  time-domain  Green’s  function,  K{r,z,t)  the 
time-domain  sensitivity  kernel  and  6c{r,  z)  the  sound  speed  perturbation. 

A  single  arrival  (arrival  G  in  Fig.  3.10)  is  considered.  The  arrival  is  refracted 
and  bottom-reflected.  A  sound  speed  perturbation  uniform  in  range  and  depth  is 
introduced  so  that  the  arrival  pattern  is  nearly  uniformly  advanced  or  delayed. 
The  linearity  of  several  different  sensitivity  kernels  is  tested:  the  full-wave  kernel 
as  given  by  eq.  3.6  for  a  one-way  transmission  (unity  source  function),  the  envelope 
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kernel  for  a  one-way  transmission,  the  full-wave  kernel  for  time-reversal  with  a 
single  source,  and  the  envelope  kernel  for  time-reversal  with  a  single  source.  The 
envelope  kernel  is  given  by [25] 


^envit') 


GR{t)  KR{t)  +  Gi{t)  Kj{t) 


(3.8) 


where  Gjiit)  and  G/(t)  are  the  real  and  imaginary  parts  of  the  unperturbed  base- 
banded  Green’s  function,  and  KR{t)  and  Kj{t)  are  the  real  and  imaginary  parts 
of  the  basebanded  pressure  sensitivity  kernels. 

The  results  of  the  linearity  tests  (Fig.  3.11)  show  the  highly  quadratic  na¬ 
ture  in  which  the  magnitude  of  the  arrival  reacts  to  sound  speed  changes.  For 
a  one-way  transmission,  using  the  envelope  of  the  sensitivity  kernel  increases  the 
range  of  linearity  over  which  the  first  Born  approximation  is  valid.  The  range  of 
linearity  for  time-reversal  (full-wave)  is  less  than  the  range  over  which  the  linear 
approximation  is  valid  for  a  one-way  transmission.  Time-reversal  effectively  delays 
each  path  such  that  they  all  arrive  at  the  same  time  at  the  receiver.  In  terms 
of  sensitivity  kernels,  this  results  in  a  single  sensitivity  kernel  for  the  peak  of  the 
focal  spot  which  is  characterized  by  several  overlapping  sensitivity  kernels,  one  for 
each  path.  [14]  The  reduced  linearity  exhibited  by  the  kernels  for  time-reversal  is 
a  a  result  of  this  overlap,  since  the  amplitude  of  each  multipath  will  affect  the 
kernel  of  the  combination  at  a  full-wave  level.  This  is  similar  to  the  problem  faced 
in  tomography,  where  an  arrival  peak  made  up  of  two  different  ray  paths  is  not 
generally  useful  for  inversion.  [12]  The  time-reversal  kernel,  therefore,  cannot  be 
expected  to  reproduce  the  evolution  of  the  focal  spot  for  very  long.  The  quadratic 
nature  of  the  peak  amplitude  variation  can  be  tackled  using  the  second  Born  ap¬ 
proximation,  demonstrated  by  Piperakis  et.  ah  in  weakly  range-dependent  ocean 
environments [26],  but  is  beyond  the  scope  of  this  paper.  The  optimized  source 
functions,  described  in  the  next  section,  are  derived  and  tested  within  the  range 
of  sound  speed  perturbations  over  which  pressure  perturbations  are  linear. 

Next,  the  sensitivity  kernels  are  used  to  calculate  Green’s  functions  in  order 
to  study  how  long  one  can  expect  the  assumption  of  linearity  to  hold  for  the  FAF05 
environment.  This  is  done  by  hrst  computing  a  baseline  and  a  perturbed  Green’s 
function  using  a  forward  run  of  the  parabolic  equation  model.  The  sensitivity 
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One-way  (full  wave)  Time-reversal  (full-wave) 


Figure  3.11:  Linearity  test  of  one-way  and  time-reversal  sensitivity  kernels.  Sound 
speeds  are  varied  uniformly  over  the  waveguide,  and  the  resulting  change  in  peak 
pressure  is  measured  as  returned  by  a  forward  run  (solid),  and  the  kernel  estimate 
(dashed) . 
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Figure  3.12:  Left:  Kernel  linearity  test  using  TC  chain  derived  ssp’s  (10  minutes 
apart).  Right:  Kernel  linearity  test  using  range  dependent  IW  model.  Blue  lines 
represent  unperturbed  Green’s  functions,  black  lines  the  perturbed  forward  run, 
and  red  lines  the  kernel  estimate. 
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kernel  is  then  used,  along  with  the  known  sound  speed  perturbation,  to  estimate  a 
Green’s  function  perturbation,  which  is  then  compared  to  that  obtained  using  the 
forward  run.  The  source  is  at  a  depth  of  112.6  m  and  receiver  at  a  depth  of  110 
m.  Two  such  tests  are  conducted.  The  first  uses  range-independent  sound  speed 
prohles  computed  using  the  thermistor  chain  data.  The  second  test  uses  range- 
dependent  sound  speed  prohles  obtained  from  the  internal  wave  model.  For  each 
test,  results  are  displayed  in  three  ways.  The  hrst  simply  displays  the  baseline, 
perturbed  and  kernel-estimated  Green’s  functions  (called  the  full-wave  Green’s 
function,  due  to  the  presence  of  carrier-frequency  modulation).  The  second  displays 
the  envelope  of  the  baseline  and  perturbed  Green’s  function  together  with  the 
kernel  estimate  computed  using  the  envelope  kernel  (eq.  3.8).  The  third  displays 
the  perturbation  to  the  envelope  compared  to  the  baseline  Green’s  function. 

For  the  left  panel  of  Fig.  3.12  sound  speed  prohles  are  computed  from  the 
thermistor  chain  data  located  close  to  the  VRA  (TG2),  at  roughly  the  same  time 
that  the  GTD  shown  in  Fig.  3.2  was  taken.  The  sound  speed  prohles  are  assumed 
to  be  range-independent.  The  Green’s  function  as  predicted  by  the  full- wave  kernel 
is  considerably  higher  in  amplitude  than  both  the  original  Green’s  function  and 
the  actual  perturbed  Green’s  function.  The  envelope  kernel  does  a  much  better 
job  of  predicting  the  perturbed  Green’s  function,  as  seen  in  the  middle  panel.  This 
is  in  agreement  with  the  lower  left  hgure  in  Fig.  3.11  where  the  envelope  kernel  is 
seen  to  have  a  larger  range  of  linearity.  The  lowest  panel  in  Fig.  3.12  shows  the 
perturbation  to  the  envelope  of  the  Green’s  function.  While  the  envelope  kernel  is 
able  to  predict  the  envelope  perturbation  at  a  10  minute  interval,  it  was  seen  to 
degrade  rapidly  for  longer  intervals. 

The  linearity  of  the  sensitivity  kernels  are  then  tested  along  with  the  range- 
dependent  sound  speed  prohles  taken  from  the  internal  wave  model.  The  results 
shown  in  the  right  hand  panels  of  Fig.  3.12.  The  baseline  Green’s  function  is 
the  same  as  in  Figs.  3.6,  3.10.  Again,  while  the  full-wave  kernel  quickly  loses  its 
validity,  the  use  of  the  envelope  kernel  is  able  to  extend  the  linear  regime  over 
which  perturbations  can  be  predicted. 
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3.5  Optimized  source  functions 

The  mathematical  framework  that  expresses  the  pressure  sensitivity  as  a 
function  of  the  source  transmission,  summarized  in  Sec.  IV,  allows  for  the  deriva¬ 
tion  of  optimized  source  functions  that  can  give  an  inverse  hlter-like  focal  spot, 
whose  quality  can  be  traded-off  for  robustness  to  sound  speed  fluctuations.  While 
the  complete  derivation  can  be  found  in  Raghukumar  et.  ah  [14],  a  brief  summary 
is  given  below.  The  main  focus  of  this  section  is  on  comparing  the  focal  spot  ob¬ 
tained  using  the  optimized  source  functions  with  experimental  time-reversal  data 
gathered  during  the  FAF05  experiment. 

3.5.1  Theory 

The  optimized  source  functions  are  derived  through  a  joint  optimization 
process  that  seeks  to  minimize  two  cost  functions:  the  first  being  a  minimization 
of  pressure  perturbations  at  the  receiver  due  to  the  changing  sound  speed,  with 
the  second  being  to  obtain  a  focal  spot  that  is  a  delta-function  in  depth  and  time. 

In  order  to  be  able  to  perform  the  derivation  using  linear  algebra,  the  var¬ 
ious  quantities  involved,  such  as  sensitivity  kernels,  sound  speed  prohles,  and 
source  functions,  are  discretized.  The  frequency-domain  sensitivity  kernel,  orig¬ 
inally  three-dimensional  (W  x  W  x  Ng),  where  W  =  60  is  the  number  of  depth 
points, iVj.  =  204  is  the  number  of  range  points,  and  Ng  =  10  the  number  of 
sources,  is  stacked  in  range  to  give  the  two-dimensional  discrete  sensitivity  kernel 
K(ci;)  with  dimensions  N^N^.  x  Ng.  Similarly  the  two-dimensional  sound  speed 
prohle  is  also  stacked  in  range  to  give  the  NzN^  x  1  vector  Ac.  The  pressure 
perturbation  (eq.  3.4)  is  then  written  as 

AP{uj)  =  [K(o;)  s{uj)f  Ac  (3.9) 

where  the  superscript  (■)^  denotes  matrix  transpose.  This  is  the  pressure  pertur- 
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bation  that  is  to  be  minimized  via  the  cost  function 

JiM  =  <  |AP|2  > 

=  s(ci;)'^  <  AcAc^  >  K(ci;)  s(ci;) 

=  s(o;)^  K(o;)^  C^c  K(o;)  s(o;) 

=  siuf  Cmdc  siu)  (3.10) 

where  Cdc  is  the  sound  speed  perturbation  covariance  matrix  and  the  superscript 
{■)^  denotes  matrix  complex  conjugate.  Cmdc  is  the  Ng  x  Ng  covariance  matrix 
of  the  Born-approximation  perturbed  Green’s  functions  at  frequency  u  between 
an  array  of  sources  and  a  single  receiver.  So  element  ij  of  Cmdc  is  the  covariance 
between  the  pressure  perturbation  at  the  receiver  from  a  transmission  from  source 
i  and  the  pressure  perturbation  at  the  receiver  from  a  transmission  from  source  j. 

The  second  condition  that  is  addressed  is  that  a  focal  spot  be  present.  In 
the  process  of  achieving  this,  the  inverse  hlter  was  re-derived.  While  time-reversal 
maximizes  the  acoustic  intensity  at  the  focal  spot,  the  inverse  hlter  adds  additional 
constraints  such  as  the  minimization  of  sidelobes  in  depth,  range  and  time.  The 
derivation  starts  with  a  Bayesian  linear  model 

u(a;)  =  G(a;)s(a;) -|- r  (3-11) 

where  u(a;)  is  an  objective  function,  whose  frequency-domain  representation  is 
a  delta  function  in  depth  and  time,  dehned  at  Np  depths.  G  is  the  Np  x  Ng 
matrix  of  Green’s  functions,  and  r  a  misht  vector.  Solving  eq.  3.11  is  equivalent 
to  minimizing  the  cost  function 

J2  =  r^Cr“^r +  s^Cs“^s  (3.12) 

The  joint  minimization  of  Ji  and  J2  with  respect  to  s  results  in  the  maxi¬ 
mum  a  posteriori  estimator[27,  14] 

S(o;)  =  (^I  +  iCmdc  +  G^G)-iG^u  (3.13) 

which  is  obtained  by  assuming  that  Cr  =  al  and  Cs  =  /dl.  7  is  introduced 
through  a  scaling  parameter  for  the  sensitivity  cost  function  Ji  when  adding  the  two 
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cost  functions  together  to  perform  the  joint  optimization.  The  optimized  source 
function  (eq.  3.13)  may  be  recognized  as  a  regularized  inverse  £lter,[8]  with  the 
terms  ^I+yCmdc  being  the  physics-based  regularization  term  that  adds  robustness 
to  environmental  perturbations  to  an  inverse  hlter-like  focal  spot.  These  fnnctions 
are  different  from  the  robust  time-reversal  method  proposed  by  Kim  et.  al.[5] 
where  the  cost  of  stability  was  a  focal  spot  being  smeared  ont  in  range,  depth  and 
time. 

3.5.2  Experimental  demonstration 

This  section  brings  together  experimental  time-reversal  data,  the  internal 
wave  model,  the  parabolic  eqnation  propagation  model,  and  the  optimized  sonrce 
fnnctions. 

In  order  to  compnte  the  optimized  sonrce  fnnctions,  the  various  components 
of  eq.  3.13  need  to  be  compnted.  The  covariance  matrix  Cmdc  is  the  weighted  inner 
prodnct  of  the  sensitivity  kernel  K(ci;)  with  the  weighting  vector  Cdc,  the  covariance 
matrix  of  sonnd  speed  pertnrbations.  As  ontlined  in  Sec.  IV,  an  experimentally 
gathered  range-independent  sound  speed  profile  is  used  to  conduct  the  forward 
model  runs.  Sensitivity  kernels  are  calculated  using  the  parabolic  eqnation  model, 
between  9  sonrces  and  11  receivers  over  a  1.4  kHz  bandwidth  (2.8  kHz-4.2  kHz). 
The  internal  wave  model  is  used  to  generate  an  ensemble  of  range-dependent  sonnd 
speed  prohles  from  which  the  covariance  matrix  for  sonnd  speed  pertnrbations  is 
calculated.  The  Green’s  fnnction  matrices,  G  may  either  be  calculated  using  the 
acoustic  propagation  code,  or  if  available,  taken  directly  from  experimental  data. 
The  presence  of  source  and  receive  vertical  arrays  during  the  FAF05  experiment 
allowed  for  the  use  of  experimental  Green’s  fnnctions  in  eq.  3.13.  Althongh  there 
were  29  source  elements  and  32  receive  elements  only  9  source  elements  and  11 
receive  elements  are  used  in  the  examples.  The  9  sonrces  are  chosen  to  be  every 
third  SRA  hydrophone  starting  from  the  bottom,  and  the  11  receivers  are  chosen  to 
be  every  fourth  VRA  hyrdrophone  starting  from  the  bottom.  This  is  done  in  order 
to  observe  greater  differences  between  time-reversal  and  the  inverse  hlter/optimized 
sonrce  fnnctions.  The  relative  gains  of  using  an  optimized  sonrce  transmission 
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scheme  are  diminished  with  an  increasing  number  of  source  elements  since  the 
quality  of  the  time-reversal  focal  spot  improves  signihcantly. 

The  next  quantity  that  needs  to  be  chosen  in  eq.  3.13  is  the  shape  of  the 
focal  spot  u(a;).  Initially,  u{z,t)  is  set  to  be  a  delta  function  in  the  baseband 
sampled  at  twice  the  bandwidth  (1.4  kHz).  This  is  then  centered  around  the 
carrier  frequency  (3.5  kHz)  and  Fourier-transformed  using  a  12  kHz  sampling  rate 
to  give  u(ci;),  whose  amplitude  is  adjusted  to  be  the  same  as  that  of  the  time- 
reversal  focal  spot.  For  a  given  bandwidth,  this  is  the  best  attainable  focal  spot. 
The  inverse  hlter  is  computed  using  eq.  3.13  with  0  =  7  =  0.  The  focal  spot 
produced  by  this  inverse  hlter  is  then  chosen  to  be  the  new  objective  function 
u(a;)  that  is  used  in  further  calculations.  The  use  of  this  objective  function  will 
result  in  the  focal  spot  cost  function  (eq.  3.12)  to  be  zero  for  the  inverse  hlter. 
Next,  the  scaling  parameter  a  is  chosen  such  that  the  focal  spot  cost  function 
for  time-reversal  =  1.  The  focal  spot  cost  function  for  the  optimized  source 
function  will  thus  be  conveniently  between  0  (inverse  hlter)  and  1  (time-reversal). 
Finally,  (3  and  7  are  varied  such  that  desirable  cost  function  characteristics  are 
achieved,  such  as  a  focal  spot  that  is  similar  to  the  inverse  hlter,  and  with  a  peak 
stability  greater  than  time-reversal.  Fig.  3.13  shows  a  representative  plot  of  the 
cost  functions  that  are  a  result  of  computing  optimized  source  functions  at  one  of 
the  times  a  CTD  was  taken  during  the  FAF  experiment  (133000  UTC  on  16  July 
2005).  Additional  plots  of  cost  functions  for  various  scenarios  may  be  found  in 
Raghukumar  et.  ah  [14] 

Time-reversal  focal  spots  are  computed  by  choosing  a  time-reversed  Green’s 
function  at  the  beginning  of  a  20-minute  window  as  a  source  shading  function,  and 
then  correlating  this  with  Green’s  functions  experimentally  acquired  every  minute. 
Optimized  source  transmissions  and  the  inverse  hlter  are  calculated  using  eq.  3.13, 
with  a  and  7  set  to  zero  to  produce  the  inverse  hlter.  Focal  spots  are  then  computed 
by  calculating  the  respective  source  shading  function  as  outlined  above,  and  then 
correlating  these  with  experimental  Green’s  functions  over  the  20-minute  window. 
In  the  FAF  experiment,  nine  GTD’s  were  taken  during  the  time  that  Green’s 
function  estimates  were  made.  Sensitivity  kernels  and  optimized  source  fuunctions 
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Figure  3.13:  Upper:  Sensitivity  cost  function  (Ji  in  eq.  3.10)  for  time-reversal 
(black),  inverse  filter  (green)  and  optimized  source  functions  (red),  and  Lower: 
Focal  spot  cost  function  (J2  in  eq.  3.12)  for  the  inverse  filter  and  optimized  source 
function.  The  plot  is  normalized  so  that  the  focal  spot  cost  function  for  time- 
reversal  is  hxed  at  1  and  that  of  the  inverse  hlter  is  0. 


are  calculated  using  the  range-independant  parabolic  equation  code  for  each  of 
the  CTD  prohles.  Focal  spots  are  compared  for  20-minute  windows  beginning  at 
the  time  the  CTD’s  were  taken.  This  allows  for  calculating  error  bars  for  the 
various  focal  spot  metrics  by  averaging  over  the  9  realizations.  The  linearity  of  the 
kernels  is  only  valid  for  approximately  10  minutes,  after  which  the  linear  hrst  Born 
approximation  breaks  down  and  the  sensitivity  kernels  are  no  longer  able  to  track 
the  evolution  of  the  Green’s  function  due  to  evolution  of  the  sound  speed  held. 
However,  focal  spots  are  still  compared  for  a  20-minute  window,  in  order  to  see  if 
results  change  signihcantly  as  a  result  of  the  kernels  no  longer  being  in  the  linear 
regime.  The  values  of  (3  and  7  were  typically  10^^  and  O.la  respectively.  The  large 
value  of  (3  is  chosen  due  to  a  being  on  the  order  of  10^^  due  to  the  units  used  in 
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representing  the  experimental  Green’s  functions.  These  choices  of  (3  and  7  result 
in  all  the  terms  in  eq.  3.13  that  are  inside  the  inverse  being  on  the  order  of  10^. 

Three  metrics  are  used  to  evaluate  the  quality  of  focal  spots  for  the  different 
schemes.  All  metrics  are  computed  relative  to  the  values  for  time-reversal.  The  hrst 
is  the  peak-to-sidelobe  ratio,  where  the  sidelobe  level  is  taken  to  be  the  magnitude 
of  the  highest  sidelobe  in  depth  or  time.  The  peak-to-sidelobe  ratio  is  a  measure 
of  the  sharpness  of  the  focal  spot.  The  second  metric  reflects  the  sensitivity  to 
sound  speed  perturbations,  and  measures  the  change  in  pressure  at  the  expected 
time  and  depth  of  the  focal  spot  due  to  changes  in  the  sound  speed  held.  The 
third  metric  is  the  peak  power  for  each  focusing  scheme. 


Time  (minutes) 


Figure  3.14:  Upper:  Peak-to-sidelobe  ratio.  Middle:  peak  stability  and  Lower: 
peak  power  for  the  inverse  hlter  (red  dashed)  and  optimized  source  function  (black 
dot-dash),  using  experimentally  derived  Green’s  functions.  The  curves  plot  the 
changes  with  respect  to  time-reversal. 

Fig.  3.14  shows  the  peak-to-sidelobe  (PSR)  ratio  (top  panel),  peak  stability 
(middle  panel)  and  peak  power  (lower  panel)  for  the  inverse  hlter  and  optimized 
source  function,  along  with  error  bars.  In  order  to  compare  with  time-reversal,  the 
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curves  for  time-reversal  have  been  subtracted  out  for  each  of  the  nine  20-minute 
windows,  after  which  the  mean  and  standard  error  are  calculated.  It  should  be 
pointed  out  that  dP  in  the  middle  panel  is  calculated  using  the  pressure  at  the 
time  of  the  peak  of  an  unperturbed  focal  spot.  Fluctuations  of  the  peak  will  cause 
it  to  be  shifted  in  time  as  a  result  of  which  dP  will  always  be  negative.  Hence, 
after  the  curve  for  time- reversal  has  been  subtracted  out,  a  positive  value  in  the 
middle  panel  of  Fig.  3.14  reflects  a  reduction  with  respect  to  time-reversal. 

The  optimized  source  function  has  a  PSR  that  starts  off  a  little  higher  than 
time-reversal,  and  builds  up  to  about  0.5  dB  higher  than  time-reversal  10  minutes 
into  the  window.  The  optimized  source  function  closely  follows  this  trend  until 
about  12  minutes  after  which  it  starts  to  show  the  highest  PSR.  The  peak  sensi¬ 
tivity  of  the  optimized  source  function  is  about  the  same  as  the  inverse  hlter  and 
time-reversal  until  about  5  minutes,  after  which  it  has  the  lowest  peak  sensitivity 
of  all  the  three  transmission  schemes.  The  peak  power  for  the  inverse  hlter  is  about 
7  dB  lower  than  time-reversal,  and  that  of  the  optimized  source  function  is  about 
2  dB  higher  than  the  inverse  hlter.  Time- reversal  has  the  highest  peak  power  since 
it  is  a  matched  hlter,  whose  specihc  goal  is  to  maximize  peak  power.  Fig.  3.14 
thus  shows  the  optimized  source  function  having  a  higher  peak  power  than  the 
inverse  hlter  (but  lower  than  time-reversal),  while  showing  a  reduced  sensitivity  to 
environmental  huctuations,  at  the  cost  of  a  lowered  PSR  compared  to  the  inverse 
hlter.  The  peak  power,  PSR  and  sensitivity  for  the  optimized  source  function  can 
be  controlled  by  varying  (3  and  7.  For  example,  in  order  to  increase  the  peak 
power,  or  7  can  be  reduced.  However,  this  would  come  at  the  cost  of  increased 
sensitivity  to  sound  speed  huctuations. 

While  the  diherences  seen  in  Fig.  3.14  are  small  and  the  error  bars  mostly 
overlap,  they  are  in  the  right  sense,  and  provide  insight  into  the  use  of  sensitivity 
kernels  for  optimal  signal  design.  In  order  to  verify  the  results,  the  same  exercise 
was  undertaken  using  just  the  internal  wave  model.  The  internal  wave  model  was 
used  to  generate  4  days  of  time-varying,  range-dependent  sound  speed  prohles, 
using  which  time-evolving  Green’s  functions  are  computed  for  the  hrst  10  minutes. 
Optimized  source  functions  were  then  computed,  and  compared  to  time-reversal 
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and  the  inverse  filter  over  nine  of  these  lO-minute  windows.  The  gains  obtained 
using  the  optimized  source  functions  were  similar  to  those  using  data.  It  was  found, 
however,  that  when  the  covariance  matrix  for  sound  speed  perturbations  used  to 
compute  the  optimized  source  functions  is  calculated  using  just  the  10  minutes 
of  sound  speed  prohles,  the  gains  are  much  higher.  Thus,  while  using  the  longer 
covariance  matrix  better  captures  the  internal  wave  dynamics  over  several  days, 
one  cannot  expect  this  to  result  in  sensitivity  reduction  over  the  relatively  short 
10  to  20  minute  windows  over  which  the  kernels  are  valid. 

Another  drawback  of  the  results  shown  in  Fig.  3.14  is  the  overlap  of  error 
bars  due  to  the  small  number  of  windows  available  to  compute  the  optimized  source 
functions.  This  limitation  is  due  to  the  limited  availability  of  CTD  prohles  during 
the  transmissions,  and  that  the  sensitivity  kernels  operate  in  linear  regime  for  only 
a  few  minutes  after  the  time  the  CTD  was  taken.  However,  the  results  show  that 
the  optimized  source  functions  are  able  to  produce  a  focal  spot  that  is  inverse 
hlter-like,  while  at  the  same  time  being  more  robust  to  environmentally  induced 
sound  speed  huctuations. 

3.6  Conclusions 

Pressure  sensitivity  kernels,  which  were  introduced  in  a  previous  paper [14], 
are  studied  using  experimental  data  via  a  multi-pronged  approach.  The  results 
of  an  internal  wave  model  are  compared  against  thermistor  chain  data  gathered 
during  the  FAF05  experiment.  A  parabolic  equation  code  is  then  used  to  model 
time- reversal  through  a  range-dependent  environment,  and  the  evolution  of  the 
focal  spot  is  compared  with  experimental  acoustic  data.  This  can  be  thought  of 
as  an  acoustic  test  of  the  internal  wave  model.  Sensitivity  kernels  are  then  derived 
using  environmental  data  gathered  during  the  experiment  with  the  intention  of 
predicting  the  evolution  of  the  time-reversal  focal  spot,  and  to  demonstrate  the 
use  of  optimized  source  functions.  The  hrst  goal  was  not  met  due  to  the  quadratic 
nature  by  which  the  peak  pressure  varies,  a  phenomenon  that  the  linear  sensitivity 
kernels  are  not  capable  of  capturing.  Time-reversal  only  exacerbates  the  problem 
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due  to  the  presence  of  overlapping  paths.  In  fact,  the  same  process  which  makes 
time-reversal  robust  to  sound  speed  fluctuations  also  reduces  the  linear  regime  for 
the  time-reversal  sensitivity  kernel. 

Sensitivity  kernels  were  then  used  to  derive  optimized  source  functions  that 
can  provide  a  focal  spot  snperior  to  time-reversal  by  taking  advantage  of  knowledge 
of  the  propagation  environment.  The  presence  of  adjustable  parameters  allows  for 
hne-tuning  focal  spot  properties  such  as  peak-to-sidelobe  ratio  or  robustness  to 
environmental  flnctuations.  These  optimized  source  functions  were  compared  to 
experimental  aconstic  data  gathered  dnring  a  recent  experiment,  which  validated 
their  ntility.  In  order  to  compnte  the  optimized  source  functions,  sensitivity  kernels 
are  calculated  using  a  parabolic  equation  model.  In  addition,  an  internal  wave 
model  (validated  against  environmental  data)  was  used  to  compute  the  covariance 
matrix  for  sonnd  speed  pertnrbations,  also  reqnired  to  compnte  the  optimized 
source  functions. 

The  limitations  of  the  scheme  proposed  are  several.  First,  the  internal 
wave  model  used  may  not  be  adeqnate,  in  which  case  insnfficient,  or  inaccurate 
knowledge  of  the  sonnd  speed  fluctuation  covariance  matrix  is  a  limiting  factor. 
The  same  is  trne  if  the  backgronnd  sonnd  speed  is  not  adeqnately  known,  in  which 
case  the  sensitivity  kernels  would  be  inaccurate.  If  sound  speed  pertnrbations  are 
large,  then  the  linearity  of  the  kernel  is  qnickly  violated.  Finally,  other  sources  of 
error,  such  as  surface  effects  (which  were  small  during  FAF05),  or  array  motion,  if 
unacconnted  for,  would  also  limit  the  use  of  the  optimized  source  functions. 

While  time-reversal  is  still  an  enormously  attractive  focusing  scheme  dne 
its  simplicity,  it  is  often  inadeqnate  for  applications  snch  as  nnderwater  commn- 
nications,  which  often  reqnire  eqnalization  in  order  to  mitigate  temporal  sidelobe 
levels.  On  the  other  hand,  the  use  of  optimized  source  functions  requires  the 
computation  of  sensitivity  kernels  and  the  covariance  matrix  for  sonnd  speed  per¬ 
tnrbations,  information  which  is  often  not  available.  However,  the  rich  reception 
dataset  available  in  FAF05  provided  observed  Green’s  fnnctions,  which  made  the 
use  of  these  optimized  source  functions  considerably  easier. 
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3. A  Internal  wave  modeling 

The  vertical  displacement  of  internal  waves,  r]{r,z,t)  can  be  expressed  as 
sum  of  vertical  eigenmodes 

p27r  pNmax 

V{r,z,t)  =  E  f{uji{k),  1, 9)(l)i{ui{k),  z)exp  i{kr  cos  9  —  ui{k)t)dud9 

Jo  J  LOj  ^ 

(3.14) 

where  uj  is  the  local  inertial  frequency,  N^ax  is  the  maximum  buoyancy  frequency 
and  k  and  (j)i{ui{k) ,  z)  are  the  horizontal  wave  number  and  vertical  eigenfunction 
respectively.  For  the  FAF05  experiment,  ui  at  43°  N  is  0.06  cycles  per  hour  and 
Nmax  is  31.8  cycles  per  hour.  Due  to  the  absence  of  azimuthal  information  in  the 
thermistor  chain  data,  the  vertical  displacement  as  given  by  the  model  is  integrated 
in  azimuth.  F{uj,l)  =  f{u!,  1,9) d9  is  a  zero  mean,  complex  Gaussian  random 

variable  associated  with  the  internal  wave  power  spectrum,  [5] 

\F{u,l)\^  =  bEGMp[  j  N{z)dzy H{l)B{u)  (3.15) 

where  hEcM  is  the  internal  wave  energy  parameter  (1  m),  p  is  the  water  density 
and  F[{1)  is  the  mode  amplitude  spectrum  given  by 

i/(/)  =  M(/2  + /2)-p/2 

and  M  is  a  normalization  factor  such  that  H{1)  =  1.  The  p  and  determine 

the  relative  mode  amplitude  by  power  law.  [5]  The  following  functional  spectrum 
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B{uj)  was  used,  which  is  compared  to  the  spectrum  of  detrended  thermistor  chain 
time  series  at  a  depth  of  39  m  in  Fig.  3.7. 


B{u)  =  W- 


Uj  UJ 


UJ, 


^^1/3 


'  cuj  + 


(3.17) 


where  hF  is  a  normalization  factor  such  that  B{uj)du  =  1.  The  spectrum 

decreases  with  a  slope  of  and  then  from  iVi/3  =  J  N{z)dz/3D. 

The  dispersion  relation  uji{k)  and  vertical  eigenfunction  (pi^ui^k),  z)  are  cal¬ 
culated  by  numerically  solving,  via  the  WAVE  code  [28] ,  the  differential  equation 

~N‘^{z)  —  00“^' 


^  ,2 

dz^ 


—  u. 


=  0 


(3.18) 


subject  to  the  boundary  conditions  0(0)  =  0(T*)  =  0  where  D  is  the  water  depth. 
N{z)  is  the  buoyancy  frequency,  calculated  from  a  representative  CTD  cast.  A 
range  of  50  horizontal  wavenumbers  k  is  chosen  to  represent  internal  wave  wave¬ 
lengths  ranging  from  20  m  to  1km.  For  each  wavenumber  k,  WAVE  returns  an 
eigenfrequency  uji{k)  for  mode  I  and  the  mode  shape  (f)i{uji{k),  z).  Mode  shapes  are 
calculated  for  11  modes,  in  order  to  be  consistent  with  the  number  of  temperature 
sensors  in  the  thermistor  chain.  In  order  to  compare  displacement  mode  shapes  as 
returned  by  the  model  with  those  derived  from  the  thermistor  chain  data,  the  dis¬ 
placement  modes  need  to  be  converted  to  temperature  modes.  This  conversion  is 
done  by  first  multiplying  the  displacement  in  eq.  3.14  by  dT/dz  obtained  from  the 
CTD  prohle  to  give  model  temperature  fluctuations.  Next,  the  covariance  matrix 
for  temperature  fluctuations  is  calculated  whose  eigenmodes  can  then  be  compared 
to  those  of  the  thermistor  chain  data.  Fig.  3.8  shows  a  comparison  between  the 
eigenmodes  calculated  by  the  model  and  that  computed  from  the  displacement 
time  series,  for  the  hrst  4  modes. 

The  covariance  matrix  can  be  calculated  directly  from  eq.  3.14  by  recogniz¬ 
ing  that 


r,{r,z,k)  ~  (3.19) 

i 

i 


rs_/ 


(3.20) 
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Discretizing  the  range  and  depth  quantities  gives  the  relation 

d?7(/c)  =  P  a  (3.21) 

where  P  is  the  x  /  matrix  (j)i{k,z)  and  a  is  the  I  x  1  vector  B{k)H{l). 

Nr  is  the  number  of  range  points  and  Nz  is  the  number  of  depth  points.  The 
covariance  matrix  for  temperature  fluctuations  can  then  be  calculated  as 

CdT  =<  dT  dT^  >fc  (3.22) 

where  dT  =  dr]  ■  dT/dz. 
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Chapter  4 


Temperature  Profile  Estimation 
using  Thermistor  Chain  Data 


This  chapter  demonstrates  an  iterative  inversion  scheme  that  inverts  ther¬ 
mistor  chain  data  for  a  reference  high-resolution  temperature  prohle.  Tempera¬ 
ture  fluctuations  that  are  observed  on  a  thermistor  chain  array  are  assumed  to  be 
caused  by  the  displacement  of  water  masses  from  different  depths,  whose  under¬ 
lying  temperature  prohle  is  unknown[l].  Given  a  sufficient  overlap  of  temperature 
time-series  between  successive  thermistor  chain  elements,  it  is  shown  that  the  un¬ 
derlying  temperature  prohle  can  be  extracted,  along  with  the  displacement  time- 
series  that  resulted  in  the  measured  thermistor  time-series.  Several  test  cases  are 
demonstrated,  starting  with  a  piecewise  linear  prohle  with  identical  displacements 
at  each  observation  depth.  A  measured  CTD  prohle  is  then  used,  along  with  a 
temperature  time-series  generated  using  a  linear  internal  wave  model.  Finally,  the 
algorithm  is  demonstrated  on  a  time-series  measured  during  the  course  of  making 
environmental  measurements  as  part  of  the  Focused  Acoustic  Fields  experiment  in 
2006. 

4.1  Introduction 

Temperature  huctuations  observed  on  a  thermistor  chain  array  are  a  valu¬ 
able  tool  to  measure  and  understand  an  ocean  environment.  Temperature  observa- 


89 


tions  can  be  easily  converted  to  sound  speed  estimates.  An  understanding  of  sound 
speed  fluctuations  is  of  vital  interest  in  acoustic  propagation  modeling[2],  in  order 
to  be  able  to  accurately  predict  the  acoustic  field  in  a  waveguide.  A  thermistor 
chain  deploys  relatively  easily,  and  can  provide  a  time-  and  depth-scale  for  temper¬ 
ature  observations.  In  shallow  water  environments,  their  utility  is  enhanced  due 
to  their  ability  to  span  almost  the  entire  water  column.  This  property  is  especially 
useful  in  studying  internal  waves  in  shallow  water. 

Internal  waves  cause  displacements  of  water  parcels  having  different  tem¬ 
peratures,  resulting  in  a  temperature  time-series [3,  4].  It  is  possible  to  directly 
relate  temperature  and  displacement,  given  a  time  and  depth  temperature  obser¬ 
vations,  along  with  a  reference  temperature  prohle.  Fig.  4.1  shows  an  example  of 
a  reference  temperature  prohle,  whose  displacement  results  in  a  temperature  time- 
series.  This  reference  temperature  prohle  can  be  difficult  to  measure.  Thermistor 
chain  elements  are  often  widely  spaced  (1-10  m  apart)  making  it  impossible  to 
obtain  a  high-resolution  temperature  estimate  in  depth.  Besides,  even  if  a  high- 
resolution  CTD  (conductivity-temperature-depth)  measurement  were  made  at  the 
same  location  as  the  thermistor  chain,  it  will  simply  be  one  realization  of  a  huc- 
tuating  temperature  prohle.  Thus,  there  is  a  need  for  a  method  of  estimating  the 
background  temperature  prohle,  the  repeated  displacement  of  which  results  in  the 
observations.  A  typical  temperature  prohle,  while  often  monotonic,  is  non-linear, 
resulting  in  a  non-linear  mapping  of  temperature  to  displacement,  and  vice-versa. 
Thus,  in  order  to  invert  the  temperature  time-series  for  a  reference  temperature 
prohle,  iterative  methods  are  an  attractive  option.  The  thermistor  chain,  in  such 
an  application,  can  be  seen  as  a  ‘synthetic  aperture’  array,  which  samples  a  re¬ 
peatedly  displaced  temperature  prohle. 

It  has  been  observed  by  Pinkel[5]  that  displacements  can  often  be  described 
by  a  Gaussian  distribution.  Given  a  reference  prohle,  temperatures  observed  at 
a  certain  depth  will  half  as  often  be  above  the  reference  prohle  at  that  depth,  as 
below  it.  Thus  the  median  of  temperature  observations  is  very  close  to  the  true 
reference,  as  shown  by  Pinkel  in  coastal  waters. 

This  paper  demonstrates  once  such  iterative  method,  based  on  the  Gauss- 
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Figure  4.1:  Illustration  of  the  principle  underlying  the  inversion  scheme.  Displace¬ 
ment  of  the  prohle  on  the  left-hand  panel  results  in  the  time-series  measured  in 
the  right-hand  panel 

Newton  method[6,  7],  that  performs  the  inversion.  The  method  is  described  in 
more  detail  in  the  next  section.  Having  described  the  method,  Sec  3  demonstrates 
the  iterative  inversion  algorithm  using  several  simulation  test  cases.  The  hrst  set 
of  test  cases  have  a  known  reference  temperature  prohle  that  can  be  compared 
with  the  estimate  to  judge  the  validity  of  the  method.  A  simple  test  case  is  that  of 
a  piecewise  linear  temperature  prohle  with  uniform  displacement  depths  at  several 
observation  points.  Next,  a  measured  CTD  prohle  is  used,  again  with  uniform 
displacements.  Finally,  displacements  are  modeled  using  a  linear  internal  wave 
model,  and  the  method  is  tested.  Section  4  applies  the  inversion  method  to  a 
measured  set  of  observations,  made  during  the  course  of  an  underwater  acoustics 
experiment.  Focused  Acoustic  Fields  2006  (FAF06).  Brief  details  of  the  experiment 
are  provided,  with  more  information  available  in  Song  et.  ah  [8].  While  the 
analysis  in  previous  chapters  utilized  environmental  and  acoustic  measurements 
made  during  FAF05,  the  thermistor  chain  elements  were  spaced  too  far  apart  in 
depth,  and  sampled  too  sparsely  in  time  for  the  estimation  method  used  in  this 
chapter  to  be  ehective.  Sec.  5  concludes  the  chapter  and  points  to  improvements 
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that  can  be  made. 


4.2  Method 

This  section  explains  the  iterative  techniqne  based  on  the  Gauss-Newton 
method  that  is  used  to  invert  the  temperature  observations  to  yield  an  estimate  of 
the  background  temperature  prohle. 

The  basic  principle  behind  the  inversion  scheme  is  that  the  temperature 
observations  are  caused  by  displacements  of  the  reference  prohle.  Thus,  each  ther¬ 
mistor  chain  element  measures  the  temperatures  of  water  parcels  associated  with 
water  depths  above  and  below.  If  between  two  successive  thermistor  elements,  all 
water  depths  between  the  two  are  measured  (the  time-series  is  long  enough),  then 
it  is  possible  to  unscramble  the  time-series  to  yield  an  estimate  of  the  true  prohle 
at  the  depths  spanned  by  the  thermistor  chain. 

Before  delving  into  the  method,  a  simple  observation  is  in  order  that  yields 
insight  into  the  validity  of  the  inversion.  If  there  are  Nzobs  observations  depths, 
each  with  a  Ntobs  point  long  time-series,  then  the  number  of  known  parameters 
is  Nzobs  X  Ntobs-  Let  the  unknown  true  prohle  be  determined  at  Nzt  depths. 
It  is  assumed  that  temperature  observations  are  caused  by  displacements  of  the 
water  column  at  Nzd  depths,  at  Ntd  points  in  time.  If  {N Zd  x  Ntd  +  Nzt)  < 
{Nzobs  X  Ntobs),  then  the  number  of  unknowns  is  less  than  the  number  of  knowns, 
and  the  inversion  has  a  unique  solution.  This  situation  can  be  expected  when 
variations  in  displacements  have  a  depth-  and/or  time-scale  greater  than  that  of 
temperature  variations.  For  example,  while  temperature  huctuations  might  have 
a  decorrelation  scale  of  1-5  m,  that  of  displacement  huctuations  might  be  10-20  m. 
This  assumption  is  used  in  the  test  cases  involving  a  uniform  displacement,  and 
shown  to  be  true  for  the  test  case  using  the  internal  wave  model.  The  internal  wave 
model  has  been  tuned  to  have  similar  statistics  as  the  temperature  observations  in 
FAF06,  and  hence  the  same  assumption  about  displacement  depth-scales  is  used 
in  Sec.  4. 

The  iterative  method  is  now  discussed.  Given  a  set  of  temperature  observa- 
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tions  Tobs{z,  t),  a  starting  estimate  for  the  reference  profile  T{z)  is  made.  Following 
Pinkel[5],the  starting  estimate  is  chosen  to  be  the  median  of  Tobs{z,  t),  interpolated 
at  Zt  depths.  An  estimate  of  the  displacement  held  ZobsiT)  can  then  be  made  using 
Tobs{z,  t)  and  the  T{z)  estimate.  As  mentioned  previously,  the  displacement  held  is 
assumed  to  vary  on  a  depth  scale  greater  than  Zobs,  and  hence  needs  to  be  mapped 
to  Zd  depths,  with  Nzd  <  Nzobs-  A  displacement  increment  held  dzd{z,t)  is  then 
computed  as  a  time-diherence,  along  with  a  T{z)  increment  held  dT{z)  referenced 
to  the  value  of  T{z)  at  the  bottom- most  depth.  The  vectorized  increment  held  is 
then  stored  as  the  vector  m,  ordered  as  the  displacement  increment  held,  followed 
by  the  temperature  prohle  increment  held. 


m  = 


(4.1) 


The  vector  m  can  be  used  to  resynthesize  a  temperature  time-series.  The  algo¬ 
rithm  then  consists  of  successively  incrementing  m  such  that  the  least-squares  ht 
between  the  temperature  held  produced  by  each  iteration  of  m  and  the  observed 
temperature  held  is  minimized. 

Let  the  function  that  maps  displacement  and  temperature  increments  to 
a  temperature  time-series  Test{z,t)  be  denoted  by  F(m).  Thus,  the  relation  that 
expresses  the  diherence  between  the  estimated  time-series  and  the  observations, 
AT{z,t),  and  the  increment  to  m  (m)  that  is  required  to  minimize  the  diherence 
can  be  expressed  as. 


AT 


dF 

dm 

Gm 


-m 


(4.2) 


where  G  can  be  recognized  as  the  Jacobian  of  F.  The  increment  m  can  then  be 
calculated  as 

m  =  (G^G  +  AI)-iGAT  (4.3) 

where  {■)^  is  the  conjugate  transpose,  or  Hermitian.  m  is  calculated  iteratively, 
where  for  each  iteration,  G  is  computed  by  incrementing  each  displacement  and 
T{z)  point,  and  storing  the  resulting  time-series  caused  by  an  individual  increment. 
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Thus  G  is  a  Nzobs  x  Ntobs  x  Npar  array,  where  Npar  =  N Zd  x  Ntd-  G  has  been 
stacked  to  produce  a  two-dimensional  array  prior  to  the  inversion  in  eq.  4.3.  The 
new  estimate  for  displacement  and  temperature  increments  is  then  calculated  as 

rrinew  =  m  m  (4.4) 

The  new  increment  held  now  allows  for  the  next  iteration.  In  the  test  cases 
described,  5-15  iterations  are  typically  sufficient  to  converge  to  a  solution. 

Various  metrics  and  cost  functions  are  computed  during  the  inversion  to 
gauge  the  quality  of  the  estimation  for  the  current  iteration.  These  are  described 
in  the  following  section  when  describing  the  test  cases. 

4.3  Test  case  -  Simulation 

This  section  describes  various  test  cases  in  order  to  demonstrate  the  inver¬ 
sion  procedure.  Test  cases  are  generated  by  hrst  starting  with  a  known  temperature 
prohle.  Known  displacements  are  then  generated  either  as  random  displacements 
or  using  an  internal  wave  model,  at  the  various  sensor  depths  in  the  water  column. 
These  displacements,  together  with  the  temperature  prohle,  are  used  to  generate 
the  temperature  observations  that  will  then  be  inverted.  Although  the  test  cases 
are  ordered  in  an  increasing  order  of  complexity  of  temperature  and  displacement 
prohles,  no  explicit  assumptions  are  made  during  the  inversion  that  indicate  a  prior 
knowledge  of  the  complexity,  or  lack  thereof. 

The  water  depth  in  all  cases  is  chosen  to  be  121  m,  with  11  sensors  placed 
between  13  m  to  63  m  (5  m  apart).  These  settings  are  chosen  based  on  their  being 
identical  to  those  used  in  the  FAF05  experiment. 

In  order  to  motivate  the  discussion,  all  plots  are  only  shown  for  hrst  test 
case.  For  all  the  others,  only  those  plots  which  show  unique  case-sensitive  behavior 
are  shown.  This  also  applies  to  the  next  section,  where  experimental  observations 
are  used. 
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4.3.1  Piecewise  linear  temperature  profile  with  uniform 
displacements 

The  first  test  case  deals  with  a  piecewise  linear  temperature  with  uniform 
displacements.  Displacements  at  any  depth  have  a  Gaussian  distribution,  with  a 
mean  displacement  of  0  m,  and  a  standard  deviation  of  2  m.  For  any  given  time, 
displacements  at  all  depths  are  equal.  Displacements  are  uncorrelated  in  time,  and 
sampled  every  10  minutes.  Fig.  4.2  shows  the  temperature  profile  and  displacement 
time-series  used  to  generate  the  temperature  observations  (Fig.  4.3). 

Due  to  identical  displacements  at  all  depths,  the  displacement  scale  in  depth 
is  large,  and  even  just  2  displacement  points  {Nzd  =  2)  should  be  sufficient  to 
estimate  the  profile.  The  depths  at  which  the  displacement  is  estimated  are  chosen 
to  be  5  m  and  95  m.  The  high  resolution  temperature  profile  is  estimated  at  1  m 
steps,  and  100  points  in  the  time-series  for  temperature  observations  are  considered. 
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Figure  4.2:  Known  piecewise  linear  temperature  profile  and  uniform  displacement 
time-series.  The  numbers  on  the  righthand  panel  indicate  sensor  depths.  Displace¬ 
ments  at  any  sensor  depth  are  uncorrelated  Gaussian  random  variables.  Displace¬ 
ments  at  any  give  time  are  equal  at  all  sensor  depths. 
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Figure  4.3:  Observed  temperature  profile  {°C)  for  the  test  case  with  a  piecewise 
linear  profile  and  uniform  displacements. 
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The  starting  profile  estimate  (Fig.  4.4)  is  chosen  to  be  the  temporal  mean 
of  the  observations,  which  is  then  linearly  extrapolated  to  N Zt  depths  at  which  the 
prohle  estimate  is  required. 


Figure  4.4:  First  T(z)  guess  extrapolated  to  Nzt  depths.  Circles  show  the  mean 
of  observations  at  each  depth,  and  the  blue  line  shows  the  extrapolation  to  N zt 
depths. 


This  prohle  estimate  is  then  used  along  with  the  temperature  observations 
Tobs{z,t)  to  get  an  estimate  for  the  displacement  time-series  Zobs(T,t),  shown  in 
Fig.  4.5.  This  hgure,  when  compared  to  the  righthand  panel  of  Fig.  4.2  reveals 
differences  in  the  displacement  estimate  requiring  the  use  of  the  inversion. 
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Another  measure  of  the  error  in  displacement  estimates  can  be  seen  in 
Fig.  4.6  which  is  essentially  the  mean  and  r.m.s  of  Fig.  4.5.  For  this  particular 
test  case,  the  mean  of  displacements  was  0  m,  and  the  r.m.s  2  m,  at  all  depths. 
Deviations  from  this  are  a  sign  of  an  incorrect  estimate. 


Figure  4.6:  Mean  (dashed)  and  r.m.s  (solid)  for  estimated  displacement  depths. 
For  this  test  case,  if  the  estimates  were  perfect,  the  mean  should  be  0,  and  r.m.s 
should  be  2. 

Fig.  4.7  shows  the  depthwise  correlation,  again  obtained  using  the  displace¬ 
ment  estimate.  This  hgure  is  yet  another  indication  of  the  current  error  in  estima¬ 
tion.  Due  to  the  uniform  nature  of  the  synthetic  displacements,  the  figure  should 
have  indicated  a  full  correlation  at  all  depths.  However,  this  not  so,  especially 
at  the  fringes,  where  observations  are  not  available.  As  a  result,  prohle  estimates 
cannot  be  expected  at  these  depths. 

Of  course,  displacement  time-series  are  only  available  for  synthetic  test 
cases.  In  more  realistic  cases,  an  idea  of  the  errors  in  displacement  estimates  can 
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Figure  4.7:  Depth-depth  correlation  coefficient  for  displacement  estimate.  For  this 
test  case,  if  the  estimate  were  perfect,  it  would  be  1  everywhere. 

only  be  obtained  by  once  again  mapping  the  displacement  estimate  to  a  temper¬ 
ature  time-series  (using  the  initial  prohle  estimate),  and  comparing  that  with  the 
observations.  Such  a  comparison  is  shown  in  Fig.  4.8  where  the  observations  have 
been  subtracted  out  to  show  the  error  in  the  estimate. 
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Figure  4.8:  Error  {°C)  in  estimate  of  temperature  observations.  Displacements  are 
computed  using  the  T(z)  estimate,  and  these  displacements  are  used  to  compute 
the  estimated  time  series.  The  error  between  the  estimate  and  observations  is 
plotted.  If  all  zero,  then  the  estimate  is  perfect. 
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Figs.  4.9  and  4.11  show  the  prior  and  posterior  parameter  uncertainities 
for  displacement  and  temperature  prohle  increments.  The  priors  are  hxed  before 
hand,  and  for  this  case  are  1  m  for  the  displacements  and  scales  with  depth  for  the 
temperature.  These  are  the  parameters  whose  increments  are  determined  after  the 
inversion.  If  a  solution  is  being  converged  upon,  these  covariances  get  smaller  with 
each  iteration,  as  seen  in  Fig.  4.11.  Once  again,  higher  values  below  70  m  reflect 
a  lack  of  observations  at  those  depths,  due  to  which  a  correct  estimate  will  never 
be  arrived  at  below  70  m. 
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Temperature  {°C) 


Figure  4.9:  Prior  displacement  (upper  panel)  and  temperature  increment  (lower 
panel)  rms.  This  is  hrst  set  by  the  user,  and  has  the  form  shown. 

The  next  hgure  (Fig.  4.10)  for  this  test  case  shows  the  prohle  estimate 
compared  to  the  true  (known)  prohle  and  the  starting  estimate.  Also  shown  is 
a  plot  of  the  error  with  respect  to  the  true  prohle.  It  is  seen  that  the  inversion 
works  well  for  depths  at  which  observations  are  available.  Beyond  these  depths,  the 
estimate  is  simply  extrapolated,  which  works  at  the  surface,  but  not  the  bottom. 
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A  gauge  of  the  quality  of  the  estimate  may  be  seen  in  Fig.  4.11  where  the  lack 
of  observations  below  63  m  result  in  poor  estimation  quality.  However,  the  hnal 
estimate  is  seen  to  be  accurate  until  70  m  which  marks  the  lowest  depth  from 
which  water  parcels  are  being  advected  to  the  sensor  array.  This  observation  also 
points  to  the  ‘synthetic  aperture’  -like  nature  of  the  estimation. 


Temperature  (°C) 


Figure  4.10:  Prohle  estimates  and  estimation  error  for  the  test  case  with  a  piecewise 
linear  prohle  and  uniform  displacements.  The  left-hand  panel  shows  the  starting 
estimate  (dashed),  hnal  estimate  (circles)  and  the  true  prohle  (solid).  The  right- 
hand  panel  shows  the  estimation  error  with  respect  to  the  true  prohle 
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Figure  4.11:  Posterior  displacement  increment  rms  (upper  panel)  and  temperature 
increment  rms  (lower  panel).  Low  values  show  a  gradual  convergence.  If  it  were  0, 
then  no  further  increments  need  be  made.  High  values  at  the  fringes  indicate  the 
absence  of  observations  at  those  depths. 
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4.3.2  CTD  profile  with  uniform  displacements 

This  test  case  adds  an  additional  level  of  complexity  to  the  previons  case 
by  having  a  more  realistic  temperatnre  profile,  while  still  having  nniform  displace¬ 
ments,  that  need  only  be  estimated  at  two  depths  (5  m  and  95  m).  The  temperatnre 
profile  shown  in  Fig.  4.12  has  been  obtained  from  a  CTD  cast  made  dnring  the 
conrse  of  the  FAF05  experiment.  Temperatnre  observations  (Fig.  4.13)  are  gen¬ 
erated  via  Ganssian  random  displacements  with  a  mean  displacement  of  0  m  and 
and  standard  deviation  of  2  m.  For  a  given  time,  displacements  are  identical  at  all 
depths.  Displacements  are  nncorrelated  in  time  and  are  10  minntes  apart.  Depth 
steps  for  the  temperatnre  prohle  estimate  are  set  to  be  1  m,  and  is  the  same  as 
the  previons  test  case. 


Fignre  4.12:  Known  FAF05  CTD  temperatnre  prohle  and  nniform  displacement 
time-series.  The  nnmbers  on  the  righthand  panel  indicate  sensor  depths.  Displace¬ 
ments  at  any  sensor  depth  are  nncorrelated  Ganssian  random  variables.  Displace¬ 
ments  at  any  give  time  are  eqnal  at  all  sensor  depths. 


Not  all  plots  shown  in  the  previons  test  case  are  shown  for  this  test  case, 
especially  since  the  displacement  time-series  are  almost  the  same. 

The  hnal  temperatnre  estimate  can  be  seen  in  Fig.  4.14,  along  with  an  error 
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Figure  4.13:  Observed  temperature  profile  (“O)  for  the  test  case  with  a  FAF05 
CTD  prohle  and  uniform  displacements. 

in  the  estimate.  Once  again,  the  estimation  is  seen  to  work  well  for  this  synthetic 
case. 

Fig.  4.15  shows  the  uncertainty  variance  for  displacement  and  temperature 
increments.  A  note  should  be  made  of  the  jump  in  uncertainty  around  30  m.  This 
can  be  attributed  to  the  non-linearities  that  arise  when  the  slope  of  T{z)  changes 
with  depth.  Displacements  will  be  sensitive  to  the  to  the  estimated  temperature 
gradient  at  that  depth.  The  displacement  may  not  be  right,  in  which  case  it 
samples  the  prohle  at  the  wrong  depth,  tries  to  correct  the  temperature  at  the 
wrong  depth,  and  the  iteration  may  not  converge. 
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Figure  4.14:  Profile  estimates  and  estimation  error  for  the  test  case  with  a  FAF05 
CTD  prohle  and  uniform  displacements.  The  left-hand  panel  shows  the  starting 
estimate  (dashed),  hnal  estimate  (circles)  and  the  true  prohle  (solid).  The  right- 
hand  panel  shows  the  estimation  error  with  respect  to  the  true  prohle. 
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Figure  4.15:  Posterior  displacement  increment  rms  (upper  panel)  and  temperature 
increment  rms  (lower  panel).  Low  values  show  a  gradual  convergence.  If  it  were  0, 
then  no  further  increments  need  be  made.  High  values  at  the  fringes  indicate  the 
absence  of  observations  at  those  depths. 


108 


4.3.3  CTD  profile  with  IW-induced  displacements 

The  last  test  case  involving  synthetic  data  utilizes  temperature  observations 
that  are  obtained  via  an  internal  wave  model.  Brief  details  regarding  the  internal 
wave  model  are  provided  below,  with  further  details  available  in  Raghukumar  et. 
ah  [9], 

The  same  CTD  prohle  that  was  used  in  the  previous  test  case  is  used  to 
compute  the  buoyancy  prohle,  which  in  turn  is  used  to  solve  for  the  vertical  modes 
and  the  dispersion  relation  for  internal  waves  in  the  given  environment.  For  a 
given  mode  number  and  wavenumber,  the  dispersion  relation  gives  information  as 
to  the  which  frequencies  the  internal  waves  can  exist.  The  buoyancy  prohle  is 
also  used  to  derive  the  internal  wave  spectrum  which  follows  a  power  law  and  is 
a  shallow  water  version  of  the  Garret-Munk  spectrum.  Mode  amplitudes  are  also 
chosen  using  the  Garett-Munk  model.  All  these  ingredients  are  then  combined 
as  a  Fourier  superposition  to  give  displacement  time  series  for  each  sensor  depth. 
Temperature  observations  can  then  be  easily  derived. 

The  true  (known)  temperature  prohle  is  assumed  to  stay  constant  through¬ 
out  the  entire  observation  interval.  This  assumption  may  not  always  be  true  in 
reality,  in  which  case  any  single  estimate  made  using  the  entire  time-series  may 
not  be  accurate. 

Due  to  the  higher  level  of  complexity  involved,  the  depth  steps  for  estimat¬ 
ing  the  high  resolution  temperature  prohle  has  been  reduced  to  0.5  m,  and  150 
points  in  the  time-series  are  used  for  the  inversion.  Data  are  samples  every  10 
minutes  so  this  corresponds  to  25  hours  of  data. 

Fig.  4.16  shows  the  CTD  prohle  and  displacement  time-series  used  to  derive 
the  temperature  observations  (Fig.  4.17). 
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Figure  4.16:  Known  FAF05  CTD  temperature  profile  and  modeled  internal  wave 
displacement  time-series.  The  numbers  on  the  righthand  panel  indicate  sensor 
depths. 
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Figure  4.17:  Observed  temperature  profile  {°C),  for  the  test  case  involving  a  FAF05 
CTD  profile  with  modeled  internal  wave  displacements. 
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Fig.  4.18  shows  the  depth-depth  correlation  coefficient  using  the  initial  dis¬ 
placement  estimate.  It  should  be  recalled  that  the  initial  displacement  estimate  is 
made  by  hrst  making  an  initial  estimate  for  the  temperature  prohle,  which  for  this 
case  was  the  mean  of  the  observations  linearly  extrapolated  to  N Zt  depths.  There 
is  signihcant  difference  in  correlation  structure  compared  to  Fig.  4.7,  due  to  the 
non-uniform  nature  of  displacements  for  this  test  case.  As  a  result,  displacements 
have  to  be  estimated  at  several  depths,  and  not  just  2  depths  as  was  done  in  the 
previous  test  cases.  Nine  displacement  depths  are  chosen,  10  m  apart,  from  5  m 
to  85  m.  There  is  also  the  option  to  have  custom  displacement  depths  based  on 
the  depth  correlation  scale,  with  closer-spaced  displacement  depths  near  the  ther- 
mocline  where  the  correlation  between  depths  is  lowest.  However,  equally  spaced 
displacement  depths  provided  the  best  estimate  for  this  case. 


Depth  (m) 


Figure  4.18:  Depth-depth  correlation  coefficient  for  displacement  estimate,  for 
the  test  case  involving  a  FAF05  CTD  prohle,  with  modeled  internal  wave 
displacements. 


112 


The  next  figure  (Fig.  4.19)  shows  the  temperature  anomaly  at  the  end  of 
15  iterations.  This  is  the  error  between  the  estimated  and  actual  temperature 
observations.  Maximum  error  is  seen  around  the  thermocline,  with  the  overall 
error  level  much  higher  compared  to  that  in  Fig.  4.8. 
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Figure  4.19:  Error  in  first  temperature  estimate  {°C).  Displacements  are  computed 
using  the  first  guess  for  T(z),  and  these  displacements  are  used  to  compute  the 
estimated  time  series.  The  error  between  the  estimate  and  observations  is  plotted. 
If  all  zero,  then  the  estimate  is  perfect. 
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A  jump  in  the  uncertainty  at  35  m  can  also  be  seen  in  the  hnal  temperature 
estimate  in  Fig.  4.20  where  the  prohle  estimate  shows  the  largest  error  aronnd  35  m. 
Depths  above  13  m  and  below  65  m  are  not  considered  due  to  a  lack  of  observations 
below  these  depths,  indicating  that  no  water  parcels  from  these  depths  are  being 
advected  to  nearby  sensors. 
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Fignre  4.20:  Prohle  estimates  and  estimation  error  for  the  test  case  with  a  FAF05 
CTD  prohle  and  modeled  internal  wave  displacements.  The  left-hand  panel  shows 
the  starting  estimate  (dashed),  hnal  estimate  (circles)  and  the  trne  prohle  (solid). 
The  right-hand  panel  shows  the  estimation  error  with  respect  to  the  trne  prohle. 


Fig.  4.21  shows  the  uncertainty  variance  for  displacement  and  temperature 
increments.  A  mnch  higher  level  of  nncertainty  is  seen  compared  to  Fig.  4.15,  with 
the  maximum  occuring  in  the  region  of  the  thermocline,  aronnd  35  m. 
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Figure  4.21:  Posterior  displacement  increment  rms  (upper  panel)  and  temperature 
increment  rms  (lower  panel).  Low  values  show  a  gradual  convergence.  If  it  were  0, 
then  no  further  increments  need  be  made.  High  values  at  the  fringes  indicate  the 
absence  of  observations  at  those  depths. 
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4.4  Inversion  applied  to  FAF06  observations 

This  section  demonstrates  the  ability  of  the  inversion  technique  to  estimate 
the  reference  temperature  profile  given  a  set  of  temperature  observations  gathered 
during  the  course  of  the  FAF06  experiment. 

A  16-element  thermistor  chain  was  deployed  over  a  period  of  17  days.  In 
order  to  demonstrate  the  inversion,  only  300  samples,  or  50  minutes  of  data  are 
used.  The  water  depth  was  121  m,  and  the  sampling  interval  10  seconds.  The 
array  spanned  water  depths  from  10  m  to  40  m.  Fig.  4.22  shows  the  temperature 
observations  over  the  period  of  interest. 
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Figure  4.22:  Observed  temperature  profile  {°C)  for  the  FAF06  experiment. 

Fig.  4.23  shows  the  displacement  estimate  computed  using  the  median  of 
temperature  observations  as  a  reference  profile.  The  presence  of  high-frequency 
content  in  the  displacement  time-series  is  an  indication  of  inaccuracies  in  the  esti¬ 
mate.  This  figure,  compared  to  Fig.  4.28  further  illustrates  this  point. 
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Fig.  4.18  shows  the  depth-depth  correlation  coefficient  using  the  initial  dis¬ 
placement  estimate.  The  correlation  between  depths  is  even  less  than  what  is 
observed  in  Fig.  4.18.  As  a  result,  custom  displacement  depths  are  used  where 
displacement  estimates  are  made  every  5  m  between  5  m  and  30  m,  and  then  every 
10  m  between  40  m  and  60  m. 


Depth  (m) 

Figure  4.24:  Depth-depth  correlation  coefficient  for  displacement  estimate,  for  the 
FAF06  thermistor  chain  data. 

Fig.  4.25  shows  the  temperature  prohle  estimate  after  6  iterations,  along 
with  the  initial  estimate.  The  estimation  was  seen  to  provide  the  lowest  error  after 
6  iterations,  after  which  the  error  was  seen  to  increase  again.  This  is  possibly  the 
result  of  the  estimation  process  not  quite  hnding  a  true  minimum.  If  this  were 
true,  a  global  optimization  method  such  as  the  use  of  a  genetic  algorithm  may  be 
required  as  a  next  step  [10].  While  a  much  higher  level  of  resolution  is  seen  in  the 
hnal  estimate,  this  hgure  by  itself  provides  no  information  about  the  accuracy  of 
the  estimate.  A  note  should  be  made  of  the  steps  that  can  be  seen  in  the  hnal 
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estimate.  Since  these  occur  in  well-observed  parts  of  the  profile,  a  lack  of  overlap 
in  the  time-series  cannot  be  an  issue.  Since  the  data  are  rather  close  to  the  water 
surface,  heating  of  the  surface  is  a  likely  candidate  for  alterations  in  the  prohle.  A 
process  such  as  this  is  not  currently  incorporated  in  the  estimation  model. 


Figure  4.25:  Starting  (dashed)  and  final  (solid)  estimates  for  the  FAF06 
observations. 

Fig.  4.26  shows  the  uncertainty  variance  for  displacement  and  temperature 
increments.  The  displacement  and  temperature  prohle  increments  show  that  while 
the  estimation  method  has  been  able  to  gradually  converge  to  a  solution  between 
20  m  and  50  m,  errors  remain  between  10  m  to  20  m.  This  is  also  rehected  in  the 
step-like  behavior  in  Fig.  4.25. 

The  next  hgure  shows  the  spectra  of  displacements  at  a  mean  depth  of 
25  m,  for  displacements  computed  using  the  median  of  temperature  observations 
(hrst  estimate)  and  hnal  temperature  prohle  estimate.  While  not  shown  in  the 
hgure,  the  displacement  spectra  were  observed  to  get  increasingly  narrow  band 
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Figure  4.26:  Posterior  displacement  (upper  panel)  and  temperature  increment 
(lower  panel)  rms.  Low  values  show  a  gradual  convergence.  If  it  were  0,  then 
no  further  increments  need  be  made.  High  values  at  the  fringes  indicate  the  ab¬ 
sence  of  observations  at  those  depths. 


over  successive  iterations.  This  is  another  indication  that  the  iterative  procedure 
is  producing  realistic  displacement  profiles  since  a  narrowband  spectra  reflects  a 
smooth  displacement  profile. 

The  final  figure  (Fig.  4.28)  shows  the  modeled  displacement  time-series  us¬ 
ing  the  final  estimate.  A  comparison  with  Fig.  4.23  shows  a  smoother  displacement 
time-series,  consistent  with  Fig.  4.27. 
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Figure  4.27:  Spectrum  of  displacement  estimates  made  at  a  mean  depth  of  25  m. 
The  solid  line  shows  the  spectrum  obtained  nsing  the  median  temperatnre  prohle 
of  the  observations  in  order  to  estimate  displacements.  The  solid  line  shows  the 
spectrnm  obtained  nsing  the  final  estimate  of  temperatnre  prohle  of  the  obser¬ 
vations  in  order  to  compnte  displacements.  The  harmonics  present  in  the  initial 
estimate  are  absent  in  the  hnal  estimate. s 
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4.5  Conclusions 

This  chapter  demonstrated  the  use  of  an  inversion  technique  that  can  be 
used  to  estimate  a  high-resolution  background  temperature  prohle,  given  a  set  of 
temperature  observations.  Temperature  fluctuations  are  assumed  to  be  caused 
by  displacements  that  have  a  depth  decorrelation  scale  much  larger  than  that  of 
observed  temperature  fluctuations.  This  assumption  gives  rise  to  a  well-posed 
inversion  problem  that  is  then  iteratively  solved  using  the  Gauss-Newton  method. 

Several  test  cases  are  demonstrated,  with  displacements  ranging  from  uni¬ 
form  to  those  derived  using  an  internal  wave  model.  The  internal  wave  model  is 
matched  to  a  set  of  experimental  temperature  observations,  on  which  the  inversion 
is  hnally  performed.  Similar  results  are  seen  for  the  internal  wave  model  and  the 
experimental  observations,  which  shows  promise  for  this  method  of  inversion. 

Some  of  the  limitations  of  this  method  are  that  the  temperature  prohle  is 
required  to  be  monotonic.  In  addition,  sensors  should  be  placed  close  enough  that 
there  is  an  overlap  of  displacement  depths  between  sensors.  While  the  method 
seems  to  work  with  non-linear  prohles,  severe  non-linearities  such  as  regions  of 
rapidly  changing  slopes,  still  pose  a  bit  of  a  challenge. 
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Chapter  5 


Conclusions  and  Future  Work 


In  this  dissertation,  pressure  sensitivity  kernels  were  applied  to  time-reversal 
acoustics,  in  an  effort  to  understand  and  quantify  the  observed  stability  of  the  time- 
reversal  focal  spot.  The  time-reversal  focal  spot  has  been  previously  observed  to 
produce  a  focal  spot  over  an  hour  after  the  initial  Green’s  functions  were  acquired. 
The  focal  spot  was  produced  despite  the  presence  of  internal  wave-induced  sound 
speed  fluctuations,  which  results  in  a  change  in  the  Green’s  function  between  the 
time  it  is  acquired  and  retransmitted. 

The  first  Born,  or  single  scattering  approximation  was  used  in  deriving 
an  expression  for  a  sensitivity  kernel  that  linearly  maps  changes  in  sound  speed 
to  changes  in  pressure  at  the  receiver.  These  kernels  are  then  studied  from  the 
point  of  view  of  time-reversal  and  applied  to  sound  transmission  in  free  space  and  a 
Pekeris  waveguide.  These  simple  examples  were  useful  in  providing  insight  into  the 
stability  of  the  time  reversal  focal  spot.  It  was  found  that  sensitivity  kernels  have 
a  Fresnel-like  interference  pattern  consisting  of  regions  of  positive  and  negative 
sensitivity  to  sound  speed  changes.  During  time-reversal  transmissions,  several  of 
these  kernels,  corresponding  to  different  paths,  were  found  to  overlap  due  to  the 
fact  that  time-reversal  delays  each  path  such  that  they  all  arrive  at  the  same  time. 
Overlapping  kernels,  each  having  different  patterns  of  sensitivity,  result  in  a  net 
reduction  of  sensitivity  compared  to  the  sensitivity  of  a  single  path. 

As  a  logical  next  step,  prior  knowledge  of  the  kernels  that  make  up  the 
focal  spot  was  utilized  to  derive  optimized  source  functions  that  can  provide  an 
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environmentally-adaptive  focusing  ability.  These  functions  are  regularized  versions 
of  the  inverse  hlter,  with  a  physics-based  regularization  using  sensitivity  kernels 
and  the  covariance  of  sound  speed  fluctuations.  The  optimized  source  functions 
were  seen  to  have  a  focal  spot  quality  in-between  time-reversal  and  the  inverse 
hlter,  but  a  stability  greater  than  both.  While  time-reversal  is  still  an  enormously 
viable  focusing  scheme  due  to  its  inherent  simplicity,  the  optimized  source  functions 
may  be  useful  when  environmental  data  are  available  such  as  during  shallow-water 
tomography  experiments. 

Sensitivity  kernels  and  the  optimized  source  functions  were  then  tested  us¬ 
ing  experimental  data  gathered  during  the  FAF05  experiment.  It  was  found  that 
the  pressure  sensitivity  to  sound  speed  changes  was  quadratic  in  nature,  as  a  result 
of  which  the  linear  hrst  Born  approximation  is  only  valid  for  about  10  minutes  of 
Green’s  function  evolution.  Sensitivity  kernels  were  derived  for  the  experimen¬ 
tal  environment,  and  used  to  compute  optimized  source  functions.  Sound  speed 
statistics  for  the  environment  were  calculated  using  an  internal  wave  model  whose 
characteristics  was  matched  to  thermistor  chain  measurements.  Optimized  source 
functions  were  then  compared  to  time-reversal  and  the  inverse  hlter  and  gains 
found  to  be  in  the  right  sense,  although  small,  due  to  the  limited  time-windows 
available  during  which  the  kernels  were  linear. 

The  last  part  of  the  dissertation  demonstrated  a  method  of  inverting  ther¬ 
mistor  chain  data  for  a  reference  temperature  prohle.  Knowledge  of  a  reference 
prohle  is  useful  in  acoustic  propagation  modeling,  or  in  order  to  compute  a  sensi¬ 
tivity  kernel  to  model  pressure  and  travel-time  huctuations.  Temperature  observa¬ 
tions  were  assumed  to  be  caused  by  displacements  which  have  a  larger  time-  and 
depth-scale  compared  to  temperature.  A  Gauss-Newton  method  was  then  utilized 
to  invert  for  a  reference  temperature.  The  method  was  then  demonstrated  on  sev¬ 
eral  test  cases  where  the  estimate  was  shown  to  compare  well  with  a  known  prohle. 
Finally,  the  method  was  demonstrated  on  thermistor  chain  data  available  during 
the  FAF05  experiment,  where  it  was  shown  to  also  perform  well  with  respect  to 
several  diagnostics. 
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5.1  Future  Work 

There  are  several  aspects  of  this  dissertation  that  can  be  explored  further, 
some  of  which  may  be  summarized  as  follows. 

•  The  second  Born  approximation  may  be  applied  in  deriving  pressure  sensi¬ 
tivity  kernels,  which  will  allow  them  to  track  received  pressure  even  as  it 
varies  quadratically.  This  enhanced  tracking  ability  could  possibly  be  used 
to  predict  the  behavior  of  the  time-reveral  focal  spot  in  realistic  perturba¬ 
tions  scenarios.  This  can  also  lead  to  the  optimized  source  functions  being 
demonstrated  on  signihcantly  longer  time  scales  compared  to  the  10  minutes 
currently  demonstrated. 

•  Equalizers  in  underwater  communications  currently  track  arrival  paths  in 
order  to  minimize  intersymbol  interference.  Sensitivity  kernels  could  prove 
useful  in  enhancing  equalizer  performance  by  providing  information  on  which 
paths  are  more  likely  to  fluctuate  given  a  particular  perturbation  scenario 
compared  to  others. 

•  The  inversion  method  for  a  temperature  prohle  currently  uses  for  its  first 
estimate  the  mean,  median  or  the  hrst  of  temperature  observations.  Bet¬ 
ter  Erst  estimates  for  the  temperature  profile  could  lead  to  improved  final 
estimates. 


